From 00f7593c39d4cad03fcf61b90aa626e3474446b0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=A9lestin=20Marot?= <celestin.marot@uclouvain.be>
Date: Fri, 12 Jun 2020 09:19:04 +0200
Subject: [PATCH] clmin, clmax and clscale options integrated in hxt
---
Mesh/meshGRegionHxt.cpp | 7 +-
contrib/hxt/CMakeLists.txt | 2 +-
contrib/hxt/tetMesh/CMakeLists.txt | 2 +-
contrib/hxt/tetMesh/exe/tetMesh_CLI.c | 19 ++--
contrib/hxt/tetMesh/include/hxt_tetDelaunay.h | 6 +-
contrib/hxt/tetMesh/include/hxt_tetMesh.h | 15 ++--
.../hxt/tetMesh/include/hxt_tetNodalSize.h | 60 +++++++++++++
contrib/hxt/tetMesh/src/hxt_tetDelaunay.c | 42 ++++-----
contrib/hxt/tetMesh/src/hxt_tetMesh.c | 34 ++++---
contrib/hxt/tetMesh/src/hxt_tetNodalSize.c | 45 +++++-----
contrib/hxt/tetMesh/src/hxt_tetNodalSize.h | 29 ------
contrib/hxt/tetMesh/src/hxt_tetRefine.c | 90 +++++++------------
contrib/hxt/tetMesh/src/hxt_tetRefine.h | 6 +-
13 files changed, 183 insertions(+), 174 deletions(-)
create mode 100644 contrib/hxt/tetMesh/include/hxt_tetNodalSize.h
delete mode 100644 contrib/hxt/tetMesh/src/hxt_tetNodalSize.h
diff --git a/Mesh/meshGRegionHxt.cpp b/Mesh/meshGRegionHxt.cpp
index c2a5aa4a7c..db429e6333 100644
--- a/Mesh/meshGRegionHxt.cpp
+++ b/Mesh/meshGRegionHxt.cpp
@@ -27,7 +27,6 @@
extern "C" {
#include "hxt_tools.h"
-#include "hxt_boundary_recovery.h"
#include "hxt_tetMesh.h"
}
@@ -37,7 +36,7 @@ static HXTStatus messageCallback(HXTMessage *msg)
return HXT_STATUS_OK;
}
-static HXTStatus meshSizeCallBack(double* pts, size_t numPts, void *userData)
+static HXTStatus nodalSizesCallBack(double* pts, size_t numPts, void *userData)
{
GRegion *gr = (GRegion *)userData;
@@ -404,8 +403,8 @@ static HXTStatus _meshGRegionHxt(std::vector<GRegion *> ®ions)
0, // void* userData;
CTX::instance()->mesh.optimizeThreshold // double qualityMin;
},
- { // meshSize
- meshSizeCallBack, // HXTStatus (*callback)(double*, size_t, void* userData)
+ { // nodalSize
+ nodalSizesCallBack, // HXTStatus (*callback)(double*, size_t, void* userData)
regions[0], // void* meshSizeData; // FIXME: should be dynamic!
CTX::instance()->mesh.lcMin,
CTX::instance()->mesh.lcMax,
diff --git a/contrib/hxt/CMakeLists.txt b/contrib/hxt/CMakeLists.txt
index d69f0a0cce..bb9de91499 100644
--- a/contrib/hxt/CMakeLists.txt
+++ b/contrib/hxt/CMakeLists.txt
@@ -131,7 +131,6 @@ set(SRC
hxt_tetFlag.c
hxt_tetMesh.c
hxt_tetNodalSize.c
- hxt_tetNodalSize.h
hxt_tetOpti.c
hxt_tetOptiUtils.h
hxt_tetPartition.h
@@ -151,6 +150,7 @@ set(INC
hxt_tetDelaunay.h
hxt_tetFlag.h
hxt_tetMesh.h
+ hxt_tetNodalSize.h
hxt_tetOpti.h
hxt_tetRepair.h
hxt_vertices.h
diff --git a/contrib/hxt/tetMesh/CMakeLists.txt b/contrib/hxt/tetMesh/CMakeLists.txt
index 8071103883..2544ffe24a 100644
--- a/contrib/hxt/tetMesh/CMakeLists.txt
+++ b/contrib/hxt/tetMesh/CMakeLists.txt
@@ -43,12 +43,12 @@ set(HXT_TETMESH_SRC
"${CMAKE_CURRENT_SOURCE_DIR}/src/hxt_tetRefine.h"
"${CMAKE_CURRENT_SOURCE_DIR}/src/hxt_tetQuality.h"
"${CMAKE_CURRENT_SOURCE_DIR}/src/hxt_edgeRemoval.h"
- "${CMAKE_CURRENT_SOURCE_DIR}/src/hxt_tetNodalSize.h"
"${CMAKE_CURRENT_SOURCE_DIR}/src/hxt_tetOptiUtils.h"
"${CMAKE_CURRENT_SOURCE_DIR}/src/hxt_tetPartition.h"
"${CMAKE_CURRENT_SOURCE_DIR}/include/hxt_tetDelaunay.h"
"${CMAKE_CURRENT_SOURCE_DIR}/include/hxt_tetFlag.h"
"${CMAKE_CURRENT_SOURCE_DIR}/include/hxt_tetMesh.h"
+ "${CMAKE_CURRENT_SOURCE_DIR}/include/hxt_tetNodalSize.h"
"${CMAKE_CURRENT_SOURCE_DIR}/include/hxt_tetOpti.h"
"${CMAKE_CURRENT_SOURCE_DIR}/include/hxt_tetRepair.h"
"${CMAKE_CURRENT_SOURCE_DIR}/include/hxt_vertices.h"
diff --git a/contrib/hxt/tetMesh/exe/tetMesh_CLI.c b/contrib/hxt/tetMesh/exe/tetMesh_CLI.c
index b638f631e4..29422c9415 100644
--- a/contrib/hxt/tetMesh/exe/tetMesh_CLI.c
+++ b/contrib/hxt/tetMesh/exe/tetMesh_CLI.c
@@ -10,7 +10,6 @@
#include "hxt_opt.h"
#include "hxt_tools.h"
#include "hxt_omp.h"
-#include <float.h>
int main(int argc, char** argv) {
@@ -19,7 +18,7 @@ int main(int argc, char** argv) {
const char* output = NULL;
const char* geoFile = NULL;
- HXTTetMeshOptions options = {.refine=1, .optimize=1, .verbosity=1, .quality.min=0.35, .meshSize.min=0.0, .meshSize.max=DBL_MAX, .meshSize.scaling=1.0};
+ HXTTetMeshOptions options = {.refine=1, .optimize=1, .verbosity=1, .quality.min=0.35, .nodalSizes.factor=1.0};
HXT_CHECK( hxtAddOption('t', "omp_num_threads",
"Number of threads used for parallel work.\n"
@@ -50,17 +49,17 @@ int main(int argc, char** argv) {
HXT_CHECK( hxtAddOption('N', "no-improvement", "Do not optimize the mesh quality", HXT_NO_FLAG, NULL, &options.optimize) );
- HXT_CHECK( hxtAddOption('l', "clmin", "minimum mesh size", HXT_DOUBLE, &HXT_POSITIVE_RANGE, &options.meshSize.min));
- HXT_CHECK( hxtAddOption('h', "clmax", "maximum mesh size", HXT_DOUBLE, &HXT_POSITIVE_RANGE, &options.meshSize.max));
- HXT_CHECK( hxtAddOption('s', "clscale", "scale mesh size by a factor of VAL", HXT_DOUBLE, &HXT_POSITIVE_RANGE, &options.meshSize.scaling));
- HXT_CHECK( hxtAddOption('a', "aspect-ratio-min", "The threshold on the aspect-ratio used during the optimization", HXT_DOUBLE, &HXT_0_1_RANGE, &options.quality.min));
+ HXT_CHECK( hxtAddOption('L', "clmin", "minimum mesh size", HXT_DOUBLE, &HXT_POSITIVE_RANGE, &options.nodalSizes.min));
+ HXT_CHECK( hxtAddOption('H', "clmax", "maximum mesh size", HXT_DOUBLE, &HXT_POSITIVE_RANGE, &options.nodalSizes.max));
+ HXT_CHECK( hxtAddOption('S', "clscale", "scale mesh size by a factor of VAL", HXT_DOUBLE, &HXT_POSITIVE_RANGE, &options.nodalSizes.factor));
+ HXT_CHECK( hxtAddOption('A', "aspect-ratio-min", "The threshold on the aspect-ratio used during the optimization", HXT_DOUBLE, &HXT_0_1_RANGE, &options.quality.min));
HXT_CHECK( hxtAddOption('v', "verbosity",
"Verbosity level of output messages\n"
" * NUM=0: print no information\n"
" * NUM=1: print some information\n"
" * NUM=2: print all information", HXT_INT, &HXT_0_2_RANGE, &options.verbosity));
- HXT_CHECK( hxtAddOption('p', "stat", "Print timing for each portion of the program", HXT_FLAG, NULL, &options.stat));
+ HXT_CHECK( hxtAddOption('s', "stat", "Print timing for each portion of the program", HXT_FLAG, NULL, &options.stat));
HXT_CHECK( hxtAddTrailingOption("INPUT_FILE", HXT_EXISTING_FILENAME, NULL, &input));
HXT_CHECK( hxtAddTrailingOption("OUTPUT_FILE", HXT_STRING, NULL, &output));
@@ -77,10 +76,12 @@ int main(int argc, char** argv) {
printf("\ninput: %s\noutput: %s\n"
"Number of threads: %d\nDelaunay threads: %d\nOptim. threads: %d\n"
"reproducible: %s\nverbosity level: %d\nshow statistics: %s\n"
- "refine: %s\noptimize : %s\nmin aspect ratio: %.3f\n\n",
+ "refine: %s\noptimize : %s\nmin aspect ratio: %.3f\n"
+ "min. nodalSize: %g\nmax. nodalSize: %g\nnodalSize scaling: %g\n\n",
input?input:"-", output?output:"-", options.defaultThreads<=0?omp_get_max_threads():options.defaultThreads,
options.delaunayThreads, options.improveThreads,
- options.reproducible?T:F, (int) options.verbosity, options.stat?T:F, options.refine?T:F, options.optimize?T:F, options.quality.min);
+ options.reproducible?T:F, (int) options.verbosity, options.stat?T:F, options.refine?T:F, options.optimize?T:F, options.quality.min,
+ options.nodalSizes.min, options.nodalSizes.max, options.nodalSizes.factor);
}
diff --git a/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h b/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h
index 67fc922806..dfd0d57dfb 100644
--- a/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h
+++ b/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h
@@ -13,7 +13,7 @@
extern "C" {
#endif
-#include "hxt_mesh.h"
+#include "hxt_tetNodalSize.h"
#include "hxt_vertices.h"
#include "hxt_bbox.h"
@@ -34,8 +34,8 @@ typedef struct {
* - if bbox==NULL, the bbox is recomputed internally;
* - if bbox!=NULL, bbox must contain all vertices */
- double* nodalSizes; /**<
- * - if nodalSize==NULL, doesn't restrict nodalSize;
+ HXTNodalSizes* nodalSizes; /**<
+ * - if nodalSizes==NULL, doesn't restrict nodalSize;
* - if nodalSize!=NULL, nodalSize contains the minimum
* mesh size at each vertex.\n
* If the insertion of a vertex create an edge smaller than
diff --git a/contrib/hxt/tetMesh/include/hxt_tetMesh.h b/contrib/hxt/tetMesh/include/hxt_tetMesh.h
index 06b5e90de4..afdd0cd76e 100644
--- a/contrib/hxt/tetMesh/include/hxt_tetMesh.h
+++ b/contrib/hxt/tetMesh/include/hxt_tetMesh.h
@@ -63,14 +63,17 @@ typedef struct {
void* userData; /* user pointer given to callback */
- /* mesh size s0 at a new point is first clamped between min and max,
- then scaled by scaling. If another point of the mesh with mesh size
- s1 (also clamped and scaled) is closer than 0.5*(s0+s1), then the new point
- is filtered out */
+ /* if we want to insert a point p0 with size s0, and there is a point p1
+ with size s1 at a distance d, we check if
+
+ d > fmax(nodalSize.min, fmin(nodalSize.min, 0.5*(s0+s1))) * scaling
+
+ if it is not the case, the point is filtered out
+ */
double min; /* default value: 0.0 */
double max; /* default value: DBL_MAX (converted to DBL_MAX if <=0.0) */
- double scaling; /* default value: 1.0 (converted to 1.0 if <=0.0) */
- } meshSize;
+ double factor; /* default value: 1.0 (converted to 1.0 if <=0.0) */
+ } nodalSizes;
} HXTTetMeshOptions;
diff --git a/contrib/hxt/tetMesh/include/hxt_tetNodalSize.h b/contrib/hxt/tetMesh/include/hxt_tetNodalSize.h
new file mode 100644
index 0000000000..f0a060ea42
--- /dev/null
+++ b/contrib/hxt/tetMesh/include/hxt_tetNodalSize.h
@@ -0,0 +1,60 @@
+// Hxt - Copyright (C)
+// 2016 - 2020 UCLouvain
+//
+// See the LICENSE.txt file for license information.
+//
+// Contributor(s):
+// Célestin Marot
+
+#ifndef HXT_TETNODALSIZE_H
+#define HXT_TETNODALSIZE_H
+
+#include "hxt_mesh.h"
+#include <math.h>
+#include <float.h>
+
+typedef struct {
+ double* array;
+ HXTStatus (*callback)(double *coord, size_t n, void* userData);
+ void* userData;
+ double min;
+ double max;
+ double factor;
+} HXTNodalSizes;
+
+HXTStatus hxtNodalSizesInit(HXTMesh* mesh, HXTNodalSizes* nodalSize);
+HXTStatus hxtNodalSizesDestroy(HXTNodalSizes* nodalSize);
+
+/// Compute sizes at vertices of the mesh from meshSizeFun
+HXTStatus hxtComputeNodalSizesFromFunction(HXTMesh* mesh, HXTNodalSizes* nodalSize);
+
+/// Compute sizes at vertices of the mesh from existing edges
+HXTStatus hxtComputeNodalSizesFromTrianglesAndLines(HXTMesh* mesh, HXTNodalSizes* nodalSize);
+HXTStatus hxtComputeNodalSizesFromMesh(HXTMesh* mesh, HXTNodalSizes* nodalSize);
+
+
+
+static inline double squareDist(double v0[3], double v1[3])
+{
+ double d0 = v1[0] - v0[0];
+ double d1 = v1[1] - v0[1];
+ double d2 = v1[2] - v0[2];
+ return d0*d0 + d1*d1 + d2*d2;
+}
+
+
+static inline int isTooClose(double s0, double s1, double squareDist, const HXTNodalSizes* ns)
+{
+ if(s0!=DBL_MAX && s1!=DBL_MAX) {
+ double meanSize = fmin(ns->max, fmax(ns->min, 0.5*(s0+s1))) * ns->factor;
+ if(squareDist < meanSize * meanSize) { // we won't enter this on overflow of meansize when max is too big...
+ return 1;
+ }
+ }
+
+ return 0;
+}
+
+
+
+#endif
diff --git a/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c b/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
index d4727a01ca..8a4584013e 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
@@ -13,6 +13,7 @@
#include "hxt_tetFlag.h"
#include "hxt_tetRepair.h"
#include "hxt_sort.h"
+#include "hxt_tetNodalSize.h"
/**
@@ -273,22 +274,11 @@ static inline HXTStatus checkTetrahedron(HXTVertex* vertices, HXTPartition* part
}
-static inline HXTStatus pointIsTooClose(const double* __restrict__ p1, const double* __restrict__ p2, double nodalSize){
- double d2 = (p1[0]-p2[0])*(p1[0]-p2[0])
- + (p1[1]-p2[1])*(p1[1]-p2[1])
- + (p1[2]-p2[2])*(p1[2]-p2[2]);
- if (d2 < /*(0.94*0.94) * */nodalSize*nodalSize){
- return HXT_STATUS_INTERNAL;
- }
-
- return HXT_STATUS_OK;
-}
-
/* if one edge of the cavity is shorter than the nodalSize, return HXT_STATUS_INTERNAL */
-static inline HXTStatus filterCavity (TetLocal* local, HXTMesh *mesh, const double *nodalSizes, const uint32_t vta)
+static inline HXTStatus filterCavity (TetLocal* local, HXTMesh *mesh, const HXTNodalSizes* nodalSizes, const uint32_t vta)
{
double *vtaCoord = mesh->vertices.coord + 4*vta;
- double vtaNodalSize = nodalSizes[vta];
+ double vtaNodalSize = nodalSizes->array[vta];
for (uint64_t i = 0 ; i< local->ball.num ; i++) {
for (unsigned j=0;j<3;j++) {
@@ -296,33 +286,35 @@ static inline HXTStatus filterCavity (TetLocal* local, HXTMesh *mesh, const doub
if (j!=3 || nodej != HXT_GHOST_VERTEX){
double *Xj = mesh->vertices.coord + 4*nodej;
- double otherNodalSize = nodalSizes[nodej];
+ double otherNodalSize = nodalSizes->array[nodej];
if(otherNodalSize==DBL_MAX){
otherNodalSize = vtaNodalSize;
}
- HXT_CHECK( pointIsTooClose(vtaCoord, Xj, 0.5*( vtaNodalSize + otherNodalSize)) );
+ if( isTooClose(vtaNodalSize, otherNodalSize, squareDist(vtaCoord, Xj), nodalSizes) )
+ return HXT_STATUS_INTERNAL;
}
}
}
return HXT_STATUS_OK;
}
-static inline HXTStatus filterTet(HXTMesh* mesh, const double *nodalSizes, const uint64_t curTet, const uint32_t vta){
+static inline HXTStatus filterTet(HXTMesh* mesh, const HXTNodalSizes* nodalSizes, const uint64_t curTet, const uint32_t vta){
HXTVertex* vertices = (HXTVertex*) mesh->vertices.coord;
double *vtaCoord = vertices[vta].coord;
- double vtaNodalSize = nodalSizes[vta];
+ double vtaNodalSize = nodalSizes->array[vta];
for (unsigned j=0; j<4; j++) {
uint32_t nodej = mesh->tetrahedra.node[4*curTet+j];
if (j!=3 || nodej != HXT_GHOST_VERTEX){
double* Xj = vertices[nodej].coord;
- double otherNodalSize = nodalSizes[nodej];
+ double otherNodalSize = nodalSizes->array[nodej];
if(otherNodalSize==DBL_MAX){
otherNodalSize = vtaNodalSize;
}
- HXT_CHECK( pointIsTooClose(vtaCoord, Xj, 0.5*( vtaNodalSize + otherNodalSize)) );
+ if( isTooClose(vtaNodalSize, otherNodalSize, squareDist(vtaCoord, Xj), nodalSizes) )
+ return HXT_STATUS_INTERNAL;
}
}
return HXT_STATUS_OK;
@@ -1228,7 +1220,7 @@ static inline HXTStatus fillingACavity(HXTMesh* mesh, TetLocal* local, unsigned
static HXTStatus insertion(HXT2Sync* shared2sync,
TetLocal* local,
unsigned char* verticesID,
- const double* nodalSizes,
+ const HXTNodalSizes* nodalSizes,
uint64_t* curTet,
const uint32_t vta,
int perfectlyDelaunay){
@@ -1297,7 +1289,7 @@ static HXTStatus insertion(HXT2Sync* shared2sync,
static inline uint32_t filterOnMooreCurve(HXTVertex* vertices,
HXTNodeInfo* nodeInfo,
uint32_t n,
- double* nodalSizes)
+ HXTNodalSizes* nodalSizes)
{
uint32_t mooreSkipped = 0;
#pragma omp parallel reduction(+:mooreSkipped)
@@ -1311,8 +1303,8 @@ static inline uint32_t filterOnMooreCurve(HXTVertex* vertices,
uint32_t lastNode = nodeInfo[i].node;
if(nodeInfo[i].status==HXT_STATUS_TRYAGAIN){
double* p2 = vertices[lastNode].coord;
- double p2Size = nodalSizes[lastNode];
- if(p1!=NULL && pointIsTooClose(p1, p2, 0.5*(p1Size+p2Size))!=HXT_STATUS_OK){
+ double p2Size = nodalSizes->array[lastNode];
+ if(p1!=NULL && isTooClose(p1Size, p2Size, squareDist(p1, p2), nodalSizes)){
mooreSkipped++;
nodeInfo[i].status=HXT_STATUS_FALSE;
}
@@ -1395,7 +1387,7 @@ static HXTStatus parallelDelaunay3D(HXTMesh* mesh,
}
if(!noReordering) {
- double* sizesToInsert = options->nodalSizes + options->insertionFirst;
+ double* sizesToInsert = options->nodalSizes->array + options->insertionFirst;
size_t vertSize = nToInsert*sizeof(HXTVertex);
size_t sizeSize = nToInsert*sizeof(double);
@@ -1921,7 +1913,7 @@ static HXTStatus parallelDelaunay3D(HXTMesh* mesh,
// 5th: put vertices at the right indices
for (uint32_t i=firstShifted; i<mesh->vertices.num; i++) {
if(options->nodalSizes!=NULL){
- options->nodalSizes[vertices[i].padding.index] = options->nodalSizes[i];
+ options->nodalSizes->array[vertices[i].padding.index] = options->nodalSizes->array[i];
}
vertices[vertices[i].padding.index] = vertices[i];
}
diff --git a/contrib/hxt/tetMesh/src/hxt_tetMesh.c b/contrib/hxt/tetMesh/src/hxt_tetMesh.c
index d7f741157a..bcda2f5771 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetMesh.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetMesh.c
@@ -55,11 +55,11 @@ HXTStatus hxtTetMesh(HXTMesh* mesh,
if(options->defaultThreads>0) {
omp_set_num_threads(options->defaultThreads);
}
- if(options->meshSize.max<=0.0) {
- options->meshSize.max = DBL_MAX;
+ if(options->nodalSizes.max<=0.0) {
+ options->nodalSizes.max = DBL_MAX;
}
- if(options->meshSize.scaling<=0.0) {
- options->meshSize.scaling = 1.0;
+ if(options->nodalSizes.factor<=0.0) {
+ options->nodalSizes.factor = 1.0;
}
double t[8]={0};
@@ -155,22 +155,32 @@ HXTStatus hxtTetMesh(HXTMesh* mesh,
t[4] = omp_get_wtime();
if(options->refine){
- HXT_CHECK(hxtCreateNodalSize(mesh, &delOptions.nodalSizes));
-
- if(options->meshSize.callback!=NULL) {
- if(hxtComputeNodalSizeFromFunction(mesh, delOptions.nodalSizes, options->meshSize.callback, options->meshSize.userData)!=HXT_STATUS_OK) {
+ HXTNodalSizes nodalSizes = {
+ .array = NULL,
+ .callback = options->nodalSizes.callback,
+ .userData = options->nodalSizes.userData,
+ .min = options->nodalSizes.min,
+ .max = options->nodalSizes.max,
+ .factor = options->nodalSizes.factor
+ };
+ HXT_CHECK(hxtNodalSizesInit(mesh, &nodalSizes));
+
+ if(nodalSizes.callback!=NULL) {
+ if(hxtComputeNodalSizesFromFunction(mesh, &nodalSizes)!=HXT_STATUS_OK) {
HXT_WARNING("Initial sizes are interpolated instead of being fetched from the size map");
- HXT_CHECK(hxtComputeNodalSizeFromTrianglesAndLines(mesh, delOptions.nodalSizes));
+ HXT_CHECK(hxtComputeNodalSizesFromTrianglesAndLines(mesh, &nodalSizes));
}
}
else
- HXT_CHECK(hxtComputeNodalSizeFromTrianglesAndLines(mesh, delOptions.nodalSizes));
+ HXT_CHECK(hxtComputeNodalSizesFromTrianglesAndLines(mesh, &nodalSizes));
+
+ delOptions.nodalSizes = &nodalSizes;
HXT_CHECK( setFlagsToProcessOnlyVolumesInBrep(mesh) );
- HXT_CHECK(hxtRefineTetrahedra(mesh, &delOptions, options->meshSize.callback, options->meshSize.userData));
+ HXT_CHECK( hxtRefineTetrahedra(mesh, &delOptions) );
- HXT_CHECK( hxtDestroyNodalSize(&delOptions.nodalSizes) );
+ HXT_CHECK( hxtNodalSizesDestroy(&nodalSizes) );
HXT_INFO_COND(options->verbosity>1, "Mesh refinement finished\n");
}
diff --git a/contrib/hxt/tetMesh/src/hxt_tetNodalSize.c b/contrib/hxt/tetMesh/src/hxt_tetNodalSize.c
index a2826d2764..337871915d 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetNodalSize.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetNodalSize.c
@@ -7,43 +7,42 @@
// Célestin Marot
#include "hxt_vertices.h"
+#include "hxt_tetNodalSize.h"
-HXTStatus hxtCreateNodalSize(HXTMesh* mesh, double** nodalSizes_ptr)
+
+HXTStatus hxtNodalSizesInit(HXTMesh* mesh, HXTNodalSizes* nodalSizes)
{
- HXT_CHECK(hxtAlignedMalloc(nodalSizes_ptr,mesh->vertices.num*sizeof(double)));
+ HXT_CHECK(hxtAlignedMalloc(&nodalSizes->array,mesh->vertices.num*sizeof(double)));
return HXT_STATUS_OK;
}
-HXTStatus hxtComputeNodalSizeFromFunction(HXTMesh* mesh, double* nodalSizes,
- HXTStatus (*meshSizeFun)(double* coord, size_t n,
- void* meshSizeData),
- void* meshSizeData)
+HXTStatus hxtComputeNodalSizesFromFunction(HXTMesh* mesh, HXTNodalSizes* nodalSizes)
{
#pragma omp parallel for
for (uint32_t i = 0; i<mesh->vertices.num; i++){
mesh->vertices.coord[4 * i + 3] = -DBL_MAX;
}
- HXT_CHECK( meshSizeFun(mesh->vertices.coord, mesh->vertices.num, meshSizeData) );
+ HXT_CHECK( nodalSizes->callback(mesh->vertices.coord, mesh->vertices.num, nodalSizes->userData) );
int failed = 0;
#pragma omp parallel for reduction(||:failed)
for (uint32_t i=0; i<mesh->vertices.num; i++) {
if(mesh->vertices.coord[4 * i + 3] <= 0.)
failed = 1;
- nodalSizes[i] = mesh->vertices.coord[4 * i + 3];
+ nodalSizes->array[i] = mesh->vertices.coord[4 * i + 3];
}
return failed==0 ? HXT_STATUS_OK : HXT_STATUS_FAILED;
}
-HXTStatus hxtComputeNodalSizeFromTrianglesAndLines(HXTMesh* mesh, double* nodalSizes)
+HXTStatus hxtComputeNodalSizesFromTrianglesAndLines(HXTMesh* mesh, HXTNodalSizes* nodalSizes)
{
HXTVertex* vertices = (HXTVertex*) mesh->vertices.coord;
#pragma omp parallel for
for (uint32_t i = 0; i<mesh->vertices.num; i++){
- nodalSizes[i] = 0;
+ nodalSizes->array[i] = 0;
vertices[i].padding.hilbertDist = 0; // we use that as a counter to do the average...
}
@@ -62,8 +61,8 @@ HXTStatus hxtComputeNodalSizeFromTrianglesAndLines(HXTMesh* mesh, double* nodalS
double l = sqrt ((X1[0]-X2[0])*(X1[0]-X2[0])+
(X1[1]-X2[1])*(X1[1]-X2[1])+
(X1[2]-X2[2])*(X1[2]-X2[2]));
- nodalSizes[n1] += l;
- nodalSizes[n2] += l;
+ nodalSizes->array[n1] += l;
+ nodalSizes->array[n2] += l;
}
}
}
@@ -78,32 +77,32 @@ HXTStatus hxtComputeNodalSizeFromTrianglesAndLines(HXTMesh* mesh, double* nodalS
double l = sqrt ((X1[0]-X2[0])*(X1[0]-X2[0])+
(X1[1]-X2[1])*(X1[1]-X2[1])+
(X1[2]-X2[2])*(X1[2]-X2[2]));
- nodalSizes[n1] += l;
- nodalSizes[n2] += l;
+ nodalSizes->array[n1] += l;
+ nodalSizes->array[n2] += l;
}
#pragma omp parallel for
for (uint32_t i=0; i<mesh->vertices.num; i++)
{
if(vertices[i].padding.hilbertDist == 0) {
- nodalSizes[i] = DBL_MAX;
+ nodalSizes->array[i] = DBL_MAX;
}
else {
- nodalSizes[i] /= (double) vertices[i].padding.hilbertDist;
+ nodalSizes->array[i] /= (double) vertices[i].padding.hilbertDist;
}
}
return HXT_STATUS_OK;
}
-HXTStatus hxtComputeNodalSizeFromMesh(HXTMesh* mesh, double* nodalSizes)
+HXTStatus hxtComputeNodalSizesFromMesh(HXTMesh* mesh, HXTNodalSizes* nodalSizes)
{
#pragma omp parallel for
for (uint32_t i = 0; i<mesh->vertices.num; i++){
- nodalSizes[i] = DBL_MAX;
+ nodalSizes->array[i] = DBL_MAX;
}
- // we do not take into account hereafter nodalSizes = to DBL_MAX
+ // we do not take into account hereafter nodalSizes->array = to DBL_MAX
for (uint32_t i = 0; i<mesh->tetrahedra.num; i++){
for (uint32_t j = 0; j<3; j++){
for (uint32_t k = j+1; k<4; k++){
@@ -115,8 +114,8 @@ HXTStatus hxtComputeNodalSizeFromMesh(HXTMesh* mesh, double* nodalSizes)
double l = sqrt ((X1[0]-X2[0])*(X1[0]-X2[0])+
(X1[1]-X2[1])*(X1[1]-X2[1])+
(X1[2]-X2[2])*(X1[2]-X2[2]));
- if(l<nodalSizes[n1]) nodalSizes[n1] = l;
- if(l<nodalSizes[n2]) nodalSizes[n2] = l;
+ if(l<nodalSizes->array[n1]) nodalSizes->array[n1] = l;
+ if(l<nodalSizes->array[n2]) nodalSizes->array[n2] = l;
}
}
}
@@ -124,8 +123,8 @@ HXTStatus hxtComputeNodalSizeFromMesh(HXTMesh* mesh, double* nodalSizes)
return HXT_STATUS_OK;
}
-HXTStatus hxtDestroyNodalSize(double** nodalSizes_ptr)
+HXTStatus hxtNodalSizesDestroy(HXTNodalSizes* nodalSizes)
{
- HXT_CHECK( hxtAlignedFree(nodalSizes_ptr) );
+ HXT_CHECK( hxtAlignedFree(&nodalSizes->array) );
return HXT_STATUS_OK;
}
diff --git a/contrib/hxt/tetMesh/src/hxt_tetNodalSize.h b/contrib/hxt/tetMesh/src/hxt_tetNodalSize.h
deleted file mode 100644
index e4330092ab..0000000000
--- a/contrib/hxt/tetMesh/src/hxt_tetNodalSize.h
+++ /dev/null
@@ -1,29 +0,0 @@
-// Hxt - Copyright (C)
-// 2016 - 2020 UCLouvain
-//
-// See the LICENSE.txt file for license information.
-//
-// Contributor(s):
-// Célestin Marot
-
-#ifndef HXT_TETNODALSIZE_H
-#define HXT_TETNODALSIZE_H
-
-#include "hxt_mesh.h"
-
-HXTStatus hxtCreateNodalSize(HXTMesh* mesh, double** nodalSizes_ptr);
-HXTStatus hxtDestroyNodalSize(double** nodalSize);
-
-/// Compute sizes at vertices of the mesh from meshSizeFun
-HXTStatus hxtComputeNodalSizeFromFunction(HXTMesh* mesh, double* nodalSize,
- HXTStatus (*meshSizeFun)(double *coord, size_t n,
- void* meshSizeData),
- void* meshSizeData);
-
-/// Compute sizes at vertices of the mesh from existing edges
-HXTStatus hxtComputeNodalSizeFromTrianglesAndLines(HXTMesh* mesh, double* nodalSize);
-HXTStatus hxtComputeNodalSizeFromMesh(HXTMesh* mesh, double* nodalSize);
-
-
-
-#endif
diff --git a/contrib/hxt/tetMesh/src/hxt_tetRefine.c b/contrib/hxt/tetMesh/src/hxt_tetRefine.c
index 76b9ae7999..52e7bcc539 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetRefine.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetRefine.c
@@ -10,6 +10,7 @@
#include "predicates.h"
#include "hxt_tetFlag.h"
#include "hxt_sort.h"
+#include "hxt_tetNodalSize.h"
// mark all the points which are in mesh->(points | lines | triangles)
static void markMeshPoints(HXTMesh* mesh)
@@ -78,33 +79,11 @@ HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions)
}
-static inline double square_dist(double v0[3], double v1[3])
-{
- return (v1[0] - v0[0])*(v1[0] - v0[0]) +
- (v1[1] - v0[1])*(v1[1] - v0[1]) +
- (v1[2] - v0[2])*(v1[2] - v0[2]);
-}
-
-
-static int is_too_close(double ptSize, double vtaSize, double squareDist)
-{
- if(ptSize!=DBL_MAX) {
- double meanSize = 0.5*(ptSize+vtaSize);
- if(squareDist < /*(0.94*0.94) * */meanSize * meanSize) {
- return 1;
- }
- }
-
- return 0;
-}
-
-
-
/* compute a point (center) inside the tetrahedron which is very likely to respect the nodalsize.
* return 0 if the computed point does respect the interpolated nodalsize
* return 1 if the computed point does not respect the interpolated nodalsize
* the interpolated nodalsize is placed into center[3] */
-static int getBestCenter(double p[4][4], double nodalSize[4], double center[4])
+static int getBestCenter(double p[4][4], double nodalSize[4], double center[4], HXTNodalSizes* ns)
{
double avg = 0.0;
double num = 0;
@@ -130,12 +109,12 @@ static int getBestCenter(double p[4][4], double nodalSize[4], double center[4])
double s3 = nodalSize[3]!=DBL_MAX && nodalSize[3]>0.0 ? nodalSize[3] : avg;
// (e/s)^2 (e is the norm of the edge, s is the mean nodalSize over that edge)
- double e0l2 = square_dist(p[0], p[1])/(0.25*(s0 + s1)*(s0 + s1));
- double e1l2 = square_dist(p[0], p[2])/(0.25*(s0 + s2)*(s0 + s2));
- double e2l2 = square_dist(p[0], p[3])/(0.25*(s0 + s3)*(s0 + s3));
- double e3l2 = square_dist(p[1], p[2])/(0.25*(s1 + s2)*(s1 + s2));
- double e4l2 = square_dist(p[1], p[3])/(0.25*(s1 + s3)*(s1 + s3));
- double e5l2 = square_dist(p[2], p[3])/(0.25*(s2 + s3)*(s2 + s3));
+ double e0l2 = squareDist(p[0], p[1])/(0.25*(s0 + s1)*(s0 + s1));
+ double e1l2 = squareDist(p[0], p[2])/(0.25*(s0 + s2)*(s0 + s2));
+ double e2l2 = squareDist(p[0], p[3])/(0.25*(s0 + s3)*(s0 + s3));
+ double e3l2 = squareDist(p[1], p[2])/(0.25*(s1 + s2)*(s1 + s2));
+ double e4l2 = squareDist(p[1], p[3])/(0.25*(s1 + s3)*(s1 + s3));
+ double e5l2 = squareDist(p[2], p[3])/(0.25*(s2 + s3)*(s2 + s3));
// normally, this thing should be near 1 if the tet are well refined
// it can be below 1 only on the surface, because we averaged the line lengths to get the nodalsize
@@ -199,10 +178,10 @@ static int getBestCenter(double p[4][4], double nodalSize[4], double center[4])
center[3] = bary0*s0 + bary1*s1 + bary2*s2 + bary3*s3; // the interpolated nodalSize
circumcenterTooClose = (
- is_too_close(s0, center[3], square_dist(p[0], center)) ||
- is_too_close(s1, center[3], square_dist(p[1], center)) ||
- is_too_close(s2, center[3], square_dist(p[2], center)) ||
- is_too_close(s3, center[3], square_dist(p[3], center))
+ isTooClose(s0, center[3], squareDist(p[0], center), ns) ||
+ isTooClose(s1, center[3], squareDist(p[1], center), ns) ||
+ isTooClose(s2, center[3], squareDist(p[2], center), ns) ||
+ isTooClose(s3, center[3], squareDist(p[3], center), ns)
);
}
else { // we should also enter this if a bary is not finite
@@ -244,10 +223,10 @@ static int getBestCenter(double p[4][4], double nodalSize[4], double center[4])
otherCenter[3] = bary0*s0 + bary1*s1 + bary2*s2 + bary3*s3; // the interpolated nodalSize
otherCenterTooClose = (
- is_too_close(s0, otherCenter[3], square_dist(p[0], otherCenter)) ||
- is_too_close(s1, otherCenter[3], square_dist(p[1], otherCenter)) ||
- is_too_close(s2, otherCenter[3], square_dist(p[2], otherCenter)) ||
- is_too_close(s3, otherCenter[3], square_dist(p[3], otherCenter))
+ isTooClose(s0, otherCenter[3], squareDist(p[0], otherCenter), ns) ||
+ isTooClose(s1, otherCenter[3], squareDist(p[1], otherCenter), ns) ||
+ isTooClose(s2, otherCenter[3], squareDist(p[2], otherCenter), ns) ||
+ isTooClose(s3, otherCenter[3], squareDist(p[3], otherCenter), ns)
);
if(circumcenterOutside || !otherCenterTooClose) {
@@ -264,7 +243,7 @@ static int getBestCenter(double p[4][4], double nodalSize[4], double center[4])
/* just fill tetCoord and tetNodalSize with the coordinates and nodal size of each node of
* tetrahedron tet */
-static inline void getTetCoordAndNodalSize(HXTMesh* mesh, double* nodalSizes, uint64_t tet,
+static inline void getTetCoordAndNodalSize(HXTMesh* mesh, double* nsarray, uint64_t tet,
double tetCoord[4][4], double tetNodalSize[4])
{
uint32_t* nodes = &mesh->tetrahedra.node[4*tet];
@@ -273,7 +252,7 @@ static inline void getTetCoordAndNodalSize(HXTMesh* mesh, double* nodalSizes, ui
tetCoord[i][1] = mesh->vertices.coord[4 * nodes[i] + 1];
tetCoord[i][2] = mesh->vertices.coord[4 * nodes[i] + 2];
- tetNodalSize[i] = nodalSizes[nodes[i]];
+ tetNodalSize[i] = nsarray[nodes[i]];
}
}
@@ -379,10 +358,8 @@ static HXTStatus balanceRefineWork(HXTMesh* mesh, uint32_t* startPt, size_t* sta
}
-HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
- HXTStatus (*meshSizeFun)(double* coord, size_t n,
- void* meshSizeData),
- void* meshSizeData)
+HXTStatus hxtRefineTetrahedra(HXTMesh* mesh,
+ HXTDelaunayOptions* delOptions)
{
int maxThreads = omp_get_max_threads();
@@ -419,10 +396,11 @@ HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
double p[4][4];
double s[4];
- getTetCoordAndNodalSize(mesh, delOptions->nodalSizes, tet, p, s);
+ getTetCoordAndNodalSize(mesh, delOptions->nodalSizes->array, tet, p, s);
// TODO: do a SIMD getBestCenter function if it gets slow
- if(getBestCenter(p, s, &newVertices[4*ptIndex]) && meshSizeFun==NULL)
+ if(getBestCenter(p, s, &newVertices[4*ptIndex], delOptions->nodalSizes) &&
+ delOptions->nodalSizes->callback==NULL)
newVertices[4*ptIndex+3] = -DBL_MAX;
#ifndef NDEBUG
@@ -435,9 +413,9 @@ HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
}
}
- // second step (meshSizeFun): compute the effective mesh size at all these newly create points
- if(meshSizeFun!=NULL)
- HXT_CHECK( meshSizeFun(newVertices, totNewPts, meshSizeData) );
+ // second step (meshSizeCB): compute the effective mesh size at all these newly create points
+ if(delOptions->nodalSizes->callback!=NULL)
+ HXT_CHECK( delOptions->nodalSizes->callback(newVertices, totNewPts, delOptions->nodalSizes->userData) );
HXTStatus status;
@@ -459,15 +437,15 @@ HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
uint64_t tetID = ptToTet[i];
double p[4][4];
double s[4];
- getTetCoordAndNodalSize(mesh, delOptions->nodalSizes, tetID, p, s);
+ getTetCoordAndNodalSize(mesh, delOptions->nodalSizes->array, tetID, p, s);
double* vtaCoord = &newVertices[4*i];
double vtaSize = newVertices[4*i + 3];
-
- if(is_too_close(s[0], vtaSize, square_dist(p[0], vtaCoord)) ||
- is_too_close(s[1], vtaSize, square_dist(p[1], vtaCoord)) ||
- is_too_close(s[2], vtaSize, square_dist(p[2], vtaCoord)) ||
- is_too_close(s[3], vtaSize, square_dist(p[3], vtaCoord))){
+
+ if(isTooClose(s[0], vtaSize, squareDist(p[0], vtaCoord), delOptions->nodalSizes) ||
+ isTooClose(s[1], vtaSize, squareDist(p[1], vtaCoord), delOptions->nodalSizes) ||
+ isTooClose(s[2], vtaSize, squareDist(p[2], vtaCoord), delOptions->nodalSizes) ||
+ isTooClose(s[3], vtaSize, squareDist(p[3], vtaCoord), delOptions->nodalSizes)){
newVertices[4*i + 3] = -DBL_MAX;
}
else {
@@ -496,7 +474,7 @@ HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
if(mesh->vertices.num > mesh->vertices.size){
status=hxtAlignedRealloc(&mesh->vertices.coord, sizeof(double)*4*mesh->vertices.num);
if(status==HXT_STATUS_OK){
- status=hxtAlignedRealloc(&delOptions->nodalSizes, sizeof(double)*mesh->vertices.num);
+ status=hxtAlignedRealloc(&delOptions->nodalSizes->array, sizeof(double)*mesh->vertices.num);
mesh->vertices.size = mesh->vertices.num;
}
}
@@ -512,7 +490,7 @@ HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
mesh->vertices.coord[v*4 ] = newVertices[4*i + 0];
mesh->vertices.coord[v*4+1] = newVertices[4*i + 1];
mesh->vertices.coord[v*4+2] = newVertices[4*i + 2];
- delOptions->nodalSizes[v] = newVertices[4*i + 3];
+ delOptions->nodalSizes->array[v] = newVertices[4*i + 3];
v++;
}
}
diff --git a/contrib/hxt/tetMesh/src/hxt_tetRefine.h b/contrib/hxt/tetMesh/src/hxt_tetRefine.h
index 3d2e8447c4..c70fb03650 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetRefine.h
+++ b/contrib/hxt/tetMesh/src/hxt_tetRefine.h
@@ -15,10 +15,6 @@
HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
/// Add points at tets circumcenter in order to fullfill a mesh size constraint
-HXTStatus hxtRefineTetrahedra(HXTMesh* mesh,
- HXTDelaunayOptions* delOptions,
- HXTStatus (*meshSizeFun)(double* coord, size_t n,
- void* meshSizeData),
- void* meshSizeData);
+HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
#endif
--
GitLab