From 0e519f2d21d2a6be8122e03d38c532934937a2dc Mon Sep 17 00:00:00 2001
From: jf remacle <jean-francois.remacle@uclouvain.be>
Date: Mon, 19 Nov 2018 16:08:09 +0100
Subject: [PATCH] mesh size field in hxtMesh3D
---
Mesh/meshGRegionHxt.cpp | 9 ++++-
contrib/hxt/hxt_mesh3d.c | 72 +++++++++++++++++++++++++----------
contrib/hxt/hxt_mesh3d.h | 12 +++++-
contrib/hxt/hxt_mesh3d_main.c | 22 +++++------
contrib/hxt/hxt_mesh3d_main.h | 4 +-
5 files changed, 83 insertions(+), 36 deletions(-)
diff --git a/Mesh/meshGRegionHxt.cpp b/Mesh/meshGRegionHxt.cpp
index 6409664fae..51be474549 100644
--- a/Mesh/meshGRegionHxt.cpp
+++ b/Mesh/meshGRegionHxt.cpp
@@ -15,6 +15,7 @@
#include "MTriangle.h"
#include "MLine.h"
#include "GmshMessage.h"
+#include "BackgroundMeshTools.h"
#ifdef HAVE_HXT
extern "C" {
@@ -29,6 +30,11 @@ HXTStatus hxtGmshMsgCallback(HXTMessage* msg){
return HXT_STATUS_OK;
}
+static double hxtMeshSizeGmshCallBack (double x, double y, double z, void *userData) {
+ GRegion *gr = (GRegion*)userData;
+ double lc2 = BGM_MeshSize(gr, 0, 0, x,y,z);
+}
+
static HXTStatus getAllFacesOfAllRegions (std::vector<GRegion *> ®ions, HXTMesh *m, std::vector<GFace *> &allFaces){
std::set<GFace *, GEntityLessThan> allFacesSet;
@@ -334,7 +340,8 @@ static HXTStatus _meshGRegionHxt(std::vector<GRegion *> ®ions)
refine,
optimize,
threshold,
- hxt_boundary_recovery));
+ hxt_boundary_recovery,
+ hxtMeshSizeGmshCallBack, regions[0]));
// HXT_CHECK(hxtMeshWriteGmsh(mesh, "hxt.msh"));
diff --git a/contrib/hxt/hxt_mesh3d.c b/contrib/hxt/hxt_mesh3d.c
index b1600ec122..dc446f92a2 100644
--- a/contrib/hxt/hxt_mesh3d.c
+++ b/contrib/hxt/hxt_mesh3d.c
@@ -15,6 +15,20 @@
// }
// #endif
+HXTStatus hxtCreateNodalSizeFromFunction(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
+ double (*mesh_size)(double x, double y, double z, void* userData),
+ void* userData)
+{
+ HXT_CHECK(hxtAlignedMalloc(&delOptions->nodalSizes,mesh->vertices.num*sizeof(double)));
+
+ #pragma omp parallel for
+ for (uint32_t i=0; i<mesh->vertices.num; i++) {
+ double* coord = &mesh->vertices.coord[4*i];
+ delOptions->nodalSizes[i] = mesh_size(coord[0], coord[1], coord[2], userData);
+ }
+
+ return HXT_STATUS_OK;
+}
HXTStatus hxtCreateNodalsizeFromTrianglesAndLines(HXTMesh* mesh, HXTDelaunayOptions* delOptions)
@@ -247,7 +261,9 @@ double hxtTetCircumcenter(double a[3], double b[3], double c[3], double d[3],
}
-HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptions, HXTMeshSize* sizeField, int *nbAdd, int iter)
+static HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
+ double (*mesh_size)(double x, double y, double z, void* userData),
+ void* userData , int *nbAdd, int iter)
{
double *newVertices;
uint32_t *numCreated;
@@ -280,33 +296,46 @@ HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptio
double circumradius2 = (a[0]-circumcenter[0])*(a[0]-circumcenter[0])+
(a[1]-circumcenter[1])*(a[1]-circumcenter[1])+
(a[2]-circumcenter[2])*(a[2]-circumcenter[2]);
+
+ setProcessedFlag(mesh, i); // we do not need to refine that tetrahedra anymore
+
// all new edges will have a length equal to circumradius2
double meshSize;
// HXTStatus status = hxtMeshSizeEvaluate ( sizeField, circumcenter, &meshSize);
- double SIZES[4];
- double AVG = 0;
- int NN = 0;
- for (int j=0;j<4;j++){
- SIZES[j] = delOptions->nodalSizes[mesh->tetrahedra.node[4*i+j]];
- if (SIZES[j] != DBL_MAX){
- NN++;
- AVG += SIZES[j];
- }
+ if(u <= 0 || v <= 0 || w <= 0 || 1.-u-v-w <= 0)
+ continue;
+
+ if(mesh_size!=NULL) {
+ meshSize = mesh_size(circumcenter[0], circumcenter[1], circumcenter[2], userData);
}
- if (NN != 4){
- AVG /= NN;
+ else { // we suppose delOptions->nodalSize!=NULL
+ double SIZES[4];
+ double AVG = 0;
+ int NN = 0;
for (int j=0;j<4;j++){
- if (SIZES[j] == DBL_MAX){
- // delOptions->nodalSizes[mesh->tetrahedra.node[4*i+j]] = AVG;
- SIZES[j] = AVG;
+ SIZES[j] = delOptions->nodalSizes[mesh->tetrahedra.node[4*i+j]];
+ if (SIZES[j] != DBL_MAX){
+ NN++;
+ AVG += SIZES[j];
+ }
+ }
+ if (NN != 4){
+ AVG /= NN;
+ for (int j=0;j<4;j++){
+ if (SIZES[j] == DBL_MAX){
+ // delOptions->nodalSizes[mesh->tetrahedra.node[4*i+j]] = AVG;
+ SIZES[j] = AVG;
+ }
}
}
+
+ meshSize = SIZES[0] * (1-u-v-w) + SIZES[1] * u + SIZES[2] * v + SIZES[3] * w;
}
- meshSize = SIZES[0] * (1-u-v-w) + SIZES[1] * u + SIZES[2] * v + SIZES[3] * w;
+
- if (u > 0 && v > 0 && w > 0 && 1.-u-v-w > 0 && meshSize * meshSize * .49 < circumradius2) {
+ if (meshSize * meshSize /* .49*/ < circumradius2) {
// printf("%llu %g\n",i,sqrt(circumradius2),meshSize);
newVertices[(size_t) 4*i ] = circumcenter[0];
newVertices[(size_t) 4*i+1] = circumcenter[1];
@@ -314,8 +343,6 @@ HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptio
newVertices[(size_t) 4*i+3] = meshSize;
localAdd++;
}
-
- setProcessedFlag(mesh, i); // we do not need to refine that tetrahedra anymore
}
}
@@ -375,11 +402,14 @@ HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptio
return HXT_STATUS_OK;
}
-HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions, HXTMeshSize* sizeField) {
+HXTStatus hxtRefineTetrahedra(HXTMesh* mesh,
+ HXTDelaunayOptions* delOptions,
+ double (*mesh_size)(double x, double y, double z, void* userData),
+ void* userData) {
int iter = 0;
while(iter++ < 40){
int nbAdd=0;
- HXT_CHECK(hxtRefineTetrahedraOneStep(mesh, delOptions, sizeField, &nbAdd, iter));
+ HXT_CHECK(hxtRefineTetrahedraOneStep(mesh, delOptions, mesh_size, userData, &nbAdd, iter));
// uint32_t nb;
// HXT_CHECK(hxtColorMesh(mesh,&nb));
if (nbAdd == 0) break;
diff --git a/contrib/hxt/hxt_mesh3d.h b/contrib/hxt/hxt_mesh3d.h
index 61d5104a56..23afb959df 100644
--- a/contrib/hxt/hxt_mesh3d.h
+++ b/contrib/hxt/hxt_mesh3d.h
@@ -9,12 +9,20 @@ HXTStatus hxtCreateFaceSearchStructure(HXTMesh* mesh, uint32_t **pfaces);
//// creates a mesh with all points of the surface mesh
HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
-/// Compute sizes at vertices of the mesh from existing edges of the tetrahera
+/// Compute sizes at vertices of the mesh from a mesh_size function
+HXTStatus hxtCreateNodalSizeFromFunction(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
+ double (*mesh_size)(double x, double y, double z, void* userData),
+ void* userData);
+
+/// Compute sizes at vertices of the mesh from existing edges
HXTStatus hxtCreateNodalsizeFromTrianglesAndLines(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
HXTStatus hxtCreateNodalsizeFromMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions);
HXTStatus hxtDestroyNodalsize(HXTDelaunayOptions* delOptions);
/// Add points at tets circumcenter in order to fullfill a mesh size constraint
-HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions, HXTMeshSize* meshsize);
+HXTStatus hxtRefineTetrahedra(HXTMesh* mesh,
+ HXTDelaunayOptions* delOptions,
+ double (*mesh_size)(double x, double y, double z, void* userData),
+ void* userData);
#endif
diff --git a/contrib/hxt/hxt_mesh3d_main.c b/contrib/hxt/hxt_mesh3d_main.c
index 52011b0012..a5f3df74ff 100644
--- a/contrib/hxt/hxt_mesh3d_main.c
+++ b/contrib/hxt/hxt_mesh3d_main.c
@@ -17,7 +17,9 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
int refine,
int optimize,
double qualityThreshold,
- HXTStatus (*bnd_recovery)(HXTMesh* mesh)) {
+ HXTStatus (*bnd_recovery)(HXTMesh* mesh),
+ double (*mesh_size)(double x, double y, double z, void* userData),
+ void* userData) {
if(defaulThreads>0) {
omp_set_num_threads(defaulThreads);
@@ -100,7 +102,7 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
HXT_CHECK( hxtConstrainLinesNotInTriangles(mesh, lines2TetMap, lines2TriMap) );
HXT_CHECK( hxtColorMesh(mesh, &nbColors) );
-
+
HXT_CHECK( hxtMapColorsToBrep(mesh, nbColors, tri2TetMap) );
HXT_CHECK( hxtAlignedFree(&tri2TetMap) );
@@ -111,20 +113,18 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
if(refine){
// HXT_CHECK(hxtComputeMeshSizeFromMesh(mesh, &delOptions));
- HXT_CHECK(hxtCreateNodalsizeFromTrianglesAndLines(mesh, &delOptions));
-
+ if(mesh_size==NULL)
+ HXT_CHECK(hxtCreateNodalsizeFromTrianglesAndLines(mesh, &delOptions));
+ else
+ HXT_CHECK(hxtCreateNodalSizeFromFunction(mesh, &delOptions, mesh_size, userData) );
+
if(nbColors!=mesh->brep.numVolumes) {
HXT_CHECK( setFlagsToProcessOnlyVolumesInBrep(mesh) );
}
- HXTMeshSize *meshSize = NULL;
- // HXT_CHECK(hxtMeshSizeCreate (context,&meshSize));
- // HXT_CHECK(hxtMeshSizeCompute (meshSize, bbox.min, bbox.max, mySize, NULL));
- // printf("time from empty mesh to first insertion: %f second\n", omp_get_wtime() - time);
- HXT_CHECK(hxtRefineTetrahedra(mesh, &delOptions, meshSize));
- // HXT_CHECK(hxtMeshSizeDelete (&meshSize));
+ HXT_CHECK(hxtRefineTetrahedra(mesh, &delOptions, mesh_size, userData));
+
HXT_CHECK( hxtDestroyNodalsize(&delOptions) );
- // #endif
}
t[5] = omp_get_wtime();
diff --git a/contrib/hxt/hxt_mesh3d_main.h b/contrib/hxt/hxt_mesh3d_main.h
index 320112d674..f54e9001c6 100644
--- a/contrib/hxt/hxt_mesh3d_main.h
+++ b/contrib/hxt/hxt_mesh3d_main.h
@@ -13,6 +13,8 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
int refine,
int optimize,
double qualityThreshold,
- HXTStatus (*bnd_recovery)(HXTMesh* mesh));
+ HXTStatus (*bnd_recovery)(HXTMesh* mesh),
+ double (*mesh_size)(double x, double y, double z, void* userData),
+ void* userData);
#endif
\ No newline at end of file
--
GitLab