diff --git a/Mesh/meshGRegionHxt.cpp b/Mesh/meshGRegionHxt.cpp index 7262f578689c68bb4e3115d8d57e6fa82c878252..c2a5aa4a7cbc9aa2cd665dbde121f0760cd9e14d 100644 --- a/Mesh/meshGRegionHxt.cpp +++ b/Mesh/meshGRegionHxt.cpp @@ -31,11 +31,6 @@ extern "C" { #include "hxt_tetMesh.h" } -static HXTStatus recoveryCallback(HXTMesh *mesh, void *userData) -{ - return hxt_boundary_recovery(mesh); -} - static HXTStatus messageCallback(HXTMessage *msg) { if(msg) Msg::Info("%s", msg->string); @@ -404,13 +399,18 @@ static HXTStatus _meshGRegionHxt(std::vector<GRegion *> ®ions) 1, // int stat; 1, // int refine; CTX::instance()->mesh.optimize, // int optimize; - CTX::instance()->mesh.optimizeThreshold, // double qualityMin; - 0, // double (*qualityFun) - 0, // void* qualityData; - meshSizeCallBack, // double (*meshSizeFun) - regions[0], // void* meshSizeData; // FIXME: should be dynamic! - recoveryCallback, // HXTStatus (*recoveryFun) - 0 // void* recoveryData; + { // quality + 0, // double (*callback)(.., userData) + 0, // void* userData; + CTX::instance()->mesh.optimizeThreshold // double qualityMin; + }, + { // meshSize + meshSizeCallBack, // HXTStatus (*callback)(double*, size_t, void* userData) + regions[0], // void* meshSizeData; // FIXME: should be dynamic! + CTX::instance()->mesh.lcMin, + CTX::instance()->mesh.lcMax, + CTX::instance()->mesh.lcFactor + } }; HXT_CHECK(hxtTetMesh(mesh, &options)); diff --git a/contrib/hxt/tetMesh/CMakeLists.txt b/contrib/hxt/tetMesh/CMakeLists.txt index 29bc5b077e89bfc268c661d0b3d8c46b12c6f56c..807110388367055f81e7b2cf64ad2ebe8d995cde 100644 --- a/contrib/hxt/tetMesh/CMakeLists.txt +++ b/contrib/hxt/tetMesh/CMakeLists.txt @@ -62,13 +62,14 @@ endif() add_subdirectory(../core "${CMAKE_CURRENT_BINARY_DIR}/core") add_subdirectory(../predicates "${CMAKE_CURRENT_BINARY_DIR}/predicates") +add_subdirectory(../tetBR "${CMAKE_CURRENT_BINARY_DIR}/tetBR") target_include_directories(hxt_tetMesh PUBLIC include PRIVATE src) target_link_libraries(hxt_tetMesh PUBLIC hxt_core - PRIVATE hxt_predicates m) + PRIVATE hxt_predicates hxt_tetBR m) target_compile_features(hxt_tetMesh PRIVATE c_std_99) @@ -79,16 +80,13 @@ target_compile_features(hxt_tetMesh PRIVATE c_std_99) if(HXT_TETMESH_BUILD_CLI) # tetMesh_CLI add_executable(tetMesh_CLI "${CMAKE_CURRENT_SOURCE_DIR}/exe/tetMesh_CLI.c") + target_link_libraries(tetMesh_CLI hxt_core hxt_predicates hxt_tetBR hxt_tetMesh) - add_subdirectory(../tetBR "${CMAKE_CURRENT_BINARY_DIR}/tetBR") - - target_link_libraries(tetMesh_CLI hxt_core hxt_predicates hxt_tetMesh hxt_tetBR) - - # # Delaunay_CLI + # Delaunay_CLI add_executable(Delaunay_CLI "${CMAKE_CURRENT_SOURCE_DIR}/exe/Delaunay_CLI.c") - target_link_libraries(Delaunay_CLI hxt_core hxt_predicates hxt_tetMesh) + target_link_libraries(Delaunay_CLI hxt_core hxt_predicates hxt_tetBR hxt_tetMesh) # tetOpti_CLI add_executable(tetOpti_CLI "${CMAKE_CURRENT_SOURCE_DIR}/exe/tetOpti_CLI.c") - target_link_libraries(tetOpti_CLI hxt_core hxt_predicates hxt_tetMesh) + target_link_libraries(tetOpti_CLI hxt_core hxt_predicates hxt_tetBR hxt_tetMesh) endif() \ No newline at end of file diff --git a/contrib/hxt/tetMesh/exe/tetMesh_CLI.c b/contrib/hxt/tetMesh/exe/tetMesh_CLI.c index 088197b9a2289b8d7109af2792a9833aba1cecf4..b638f631e4710b2035d7be301c4bba12fbc16017 100644 --- a/contrib/hxt/tetMesh/exe/tetMesh_CLI.c +++ b/contrib/hxt/tetMesh/exe/tetMesh_CLI.c @@ -7,25 +7,19 @@ // CĂ©lestin Marot #include "hxt_tetMesh.h" -#include "hxt_boundary_recovery.h" #include "hxt_opt.h" #include "hxt_tools.h" #include "hxt_omp.h" +#include <float.h> -static HXTStatus recoveryFun(HXTMesh* mesh, void* userData) { - HXT_UNUSED(userData); - HXT_CHECK( hxt_boundary_recovery(mesh) ); - return HXT_STATUS_OK; -} - int main(int argc, char** argv) { const char* input = NULL; const char* output = NULL; const char* geoFile = NULL; - HXTTetMeshOptions options = {.refine=1, .optimize=1, .verbosity=1, .qualityMin=0.35, .recoveryFun=recoveryFun}; + HXTTetMeshOptions options = {.refine=1, .optimize=1, .verbosity=1, .quality.min=0.35, .meshSize.min=0.0, .meshSize.max=DBL_MAX, .meshSize.scaling=1.0}; HXT_CHECK( hxtAddOption('t', "omp_num_threads", "Number of threads used for parallel work.\n" @@ -46,20 +40,27 @@ int main(int argc, char** argv) { "if given\n" " - the same vertices in the same order\n" " - the same number of threads\n" + "and if compiled with the same compiler on the same machine.\n" "WARNING: This option slows down the mesh generation\n" "except with 1 thread", HXT_FLAG, NULL, &options.reproducible)); - 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('s', "stat", "Print timing for each portion of the program", HXT_FLAG, NULL, &options.stat)); - HXT_CHECK( hxtAddOption('n', "no-refinement", "Do not refine the mesh => no vertices are added", HXT_NO_FLAG, NULL, &options.refine)); - HXT_CHECK( hxtAddOption('N', "no-improvement", "Do not optimize the mesh quality", HXT_NO_FLAG, NULL, &options.optimize) ); - HXT_CHECK( hxtAddOption('a', "aspect-ratio-min", "The threshold on the aspect-ratio used during the optimization", HXT_DOUBLE, &HXT_0_1_RANGE, &options.qualityMin)); const HXTOptionArgumentConstraints geoSuffix = {.stringSuffix=".geo"}; HXT_CHECK( hxtAddOption('g', "generate-geo", "Write a .geo file that describe the volumes from the input surface mesh", HXT_ASK_TO_ERASE_FILENAME, &geoSuffix, &geoFile)); + HXT_CHECK( hxtAddOption('n', "no-refinement", "Do not refine the mesh => no vertices are added", HXT_NO_FLAG, NULL, &options.refine)); + 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('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( hxtAddTrailingOption("INPUT_FILE", HXT_EXISTING_FILENAME, NULL, &input)); HXT_CHECK( hxtAddTrailingOption("OUTPUT_FILE", HXT_STRING, NULL, &output)); @@ -79,7 +80,7 @@ int main(int argc, char** argv) { "refine: %s\noptimize : %s\nmin aspect ratio: %.3f\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.qualityMin); + options.reproducible?T:F, (int) options.verbosity, options.stat?T:F, options.refine?T:F, options.optimize?T:F, options.quality.min); } diff --git a/contrib/hxt/tetMesh/include/hxt_tetMesh.h b/contrib/hxt/tetMesh/include/hxt_tetMesh.h index d47aac46b75ba2398bc4fe208591bd4c2b9ab5ee..06b5e90de4b9251e683a2c42e69f729336c67d99 100644 --- a/contrib/hxt/tetMesh/include/hxt_tetMesh.h +++ b/contrib/hxt/tetMesh/include/hxt_tetMesh.h @@ -27,7 +27,8 @@ typedef struct { - `delaunayThreads<0`: use OMP_NUM_PROCS */ int improveThreads; /* Same as above but for mesh improvements operations */ int reproducible; /*The output mesh will be identical with this option - if given + if the machine and compiler are identical + and if given - the same vertices in the same order - the same number of threads WARNING: This option slows down the mesh generation @@ -44,22 +45,32 @@ typedef struct { above the threshold given by `qualityMin` (Boolean option) */ - double qualityMin; /* Threshold for mesh improvement */ + struct { + /* function giving the quality of an element, or NULL to use default function*/ + double (*callback)(double* p0, double* p1, double* p2, double* p3, void* userData); - /* function giving the quality of an element, or NULL to use default function*/ - double (*qualityFun)(double* p0,double* p1, double* p2, double* p3, - void* qualityData); - void* qualityData; /* user pointer to give to qualityFun */ + void* userData; /* user pointer given to callback */ - /* function giving the desired mesh size at a position in space, - or NULL to use default function*/ - HXTStatus (*meshSizeFun)(double* coord, size_t n, - void* meshSizeData); - void* meshSizeData; /* user pointer to give to meshSizeFun */ + double min; /* Threshold for mesh improvement (Hxt tries to improve all tet + with a quality below quality.min */ + } quality; - /* function to recover missing features in a mesh */ - HXTStatus (*recoveryFun)(HXTMesh* mesh, void* recoveryData); - void* recoveryData; /* user pointer to give to recoveryFun */ + + struct { + /* function giving the desired mesh size at a position in space, + or NULL to use default function*/ + HXTStatus (*callback)(double* coord, size_t n, void* userData); + + 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 */ + 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; } HXTTetMeshOptions; @@ -73,8 +84,7 @@ typedef struct { * .stat = 1, * .refine = 1, * .optimize = 1, - * .qualityMin = 0.35, - * .recoveryFun = hxt_boundary_recovery + * .quality.min = 0.35, * } * * Everything that is not explicitely initialized is set to zero. @@ -85,8 +95,7 @@ typedef struct { * hxtTetMesh3d(mesh, &(HXTMesh3dOptions) {.refine=1, * .optimize=1, * .stat=1, - * .qualityMin=1, - * .bndRecovery=hxt_boundary_recovery }); + * .quality.min=0.35 }); * * This is even shorter than without using a structure. */ diff --git a/contrib/hxt/tetMesh/src/hxt_tetMesh.c b/contrib/hxt/tetMesh/src/hxt_tetMesh.c index 92c21e38cf6def68a717192519765d1097197dd6..d7f741157a9a35eee6b570c6c6ae298405cc736c 100644 --- a/contrib/hxt/tetMesh/src/hxt_tetMesh.c +++ b/contrib/hxt/tetMesh/src/hxt_tetMesh.c @@ -18,6 +18,8 @@ #include "hxt_tetQuality.h" #include "hxt_omp.h" +#include "hxt_boundary_recovery.h" + // void aspect_ratio_graph(HXTMesh* mesh) { // // make a count of aspect ratio... // unsigned bins[21] = {0}; @@ -53,6 +55,12 @@ 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->meshSize.scaling<=0.0) { + options->meshSize.scaling = 1.0; + } double t[8]={0}; t[0] = omp_get_wtime(); @@ -97,10 +105,6 @@ HXTStatus hxtTetMesh(HXTMesh* mesh, t[2] = omp_get_wtime(); if (nbMissingTriangles!=0 || nbMissingLines!=0){ - if(options->recoveryFun==NULL) - return HXT_ERROR_MSG(HXT_STATUS_ERROR, - "there are missing features but no boundary recovery function is given"); - if(nbMissingTriangles) HXT_INFO("Recovering %" HXTu64 " missing facet(s)", nbMissingTriangles); @@ -109,7 +113,7 @@ HXTStatus hxtTetMesh(HXTMesh* mesh, nbMissingLines); uint32_t oldNumVertices = mesh->vertices.num; - HXT_CHECK(options->recoveryFun(mesh, options->recoveryData)); + HXT_CHECK( hxt_boundary_recovery(mesh) ); if(oldNumVertices < mesh->vertices.num) { HXT_INFO("Steiner(s) point(s) were inserted"); @@ -153,8 +157,8 @@ HXTStatus hxtTetMesh(HXTMesh* mesh, if(options->refine){ HXT_CHECK(hxtCreateNodalSize(mesh, &delOptions.nodalSizes)); - if(options->meshSizeFun!=NULL) { - if(hxtComputeNodalSizeFromFunction(mesh, delOptions.nodalSizes, options->meshSizeFun, options->meshSizeData)!=HXT_STATUS_OK) { + if(options->meshSize.callback!=NULL) { + if(hxtComputeNodalSizeFromFunction(mesh, delOptions.nodalSizes, options->meshSize.callback, options->meshSize.userData)!=HXT_STATUS_OK) { HXT_WARNING("Initial sizes are interpolated instead of being fetched from the size map"); HXT_CHECK(hxtComputeNodalSizeFromTrianglesAndLines(mesh, delOptions.nodalSizes)); } @@ -164,7 +168,7 @@ HXTStatus hxtTetMesh(HXTMesh* mesh, HXT_CHECK( setFlagsToProcessOnlyVolumesInBrep(mesh) ); - HXT_CHECK(hxtRefineTetrahedra(mesh, &delOptions, options->meshSizeFun, options->meshSizeData)); + HXT_CHECK(hxtRefineTetrahedra(mesh, &delOptions, options->meshSize.callback, options->meshSize.userData)); HXT_CHECK( hxtDestroyNodalSize(&delOptions.nodalSizes) ); @@ -181,9 +185,9 @@ HXTStatus hxtTetMesh(HXTMesh* mesh, HXT_CHECK( hxtOptimizeTetrahedra(mesh, &(HXTOptimizeOptions) { .bbox = &bbox, - .qualityFun = options->qualityFun, - .qualityData = options->qualityData, - .qualityMin = options->qualityMin, + .qualityFun = options->quality.callback, + .qualityData = options->quality.userData, + .qualityMin = options->quality.min, .numThreads = options->improveThreads, .numVerticesConstrained = numVerticesConstrained, .verbosity = options->verbosity,