mesh size field in hxtMesh3D

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 *> &regions, HXTMesh *m, std::vector<GFace *> &allFaces){
   std::set<GFace *, GEntityLessThan> allFacesSet;
@@ -334,7 +340,8 @@ static HXTStatus _meshGRegionHxt(std::vector<GRegion *> &regions)
 			 refine,
                          optimize,
 			 threshold,
-			 hxt_boundary_recovery));  
+			 hxt_boundary_recovery,
+			 hxtMeshSizeGmshCallBack, regions[0]));  
 
   
   //  HXT_CHECK(hxtMeshWriteGmsh(mesh, "hxt.msh"));

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,10 +296,20 @@ 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);
 
+        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);
+        }
+        else { // we suppose delOptions->nodalSize!=NULL
           double SIZES[4];
           double AVG = 0;
           int NN = 0;
@@ -305,8 +331,11 @@ HXTStatus hxtRefineTetrahedraOneStep(HXTMesh* mesh, HXTDelaunayOptions* delOptio
           }
 
           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;

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

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);
@@ -111,20 +113,18 @@ HXTStatus hxtTetMesh3d(HXTMesh* mesh,
 
   if(refine){
     // HXT_CHECK(hxtComputeMeshSizeFromMesh(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();

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