octree

CMakeLists.txt

@@ -78,6 +78,7 @@ opt(OPENMP "Enable OpenMP" OFF)
 opt(OPTHOM "Enable high-order mesh optimization tools" ${DEFAULT})
 opt(OS_SPECIFIC_INSTALL "Enable OS-specific (e.g. app bundle) installation" OFF)
 opt(OSMESA "Enable OSMesa for offscreen rendering (experimental)" OFF)
+opt(P4EST "Enable p4est for enabling automatic mesh size firld (experimental)" OFF)
 opt(PARSER "Enable GEO file parser (required for .geo/.pos files)" ${DEFAULT})
 opt(PETSC "Enable PETSc linear solvers (required for SLEPc)" OFF)
 opt(PLUGINS "Enable post-processing plugins" ${DEFAULT})
@@ -850,6 +851,17 @@ if(ENABLE_POPPLER)
   endif()
 endif()
 
+if(ENABLE_P4EST)
+  find_library(P4EST_LIB p4est)
+  find_path(P4EST_INC "p4est.h" PATH_SUFFIXES src)
+  if(P4EST_LIB AND P4EST_INC)
+    set_config_option(HAVE_P4EST "p4est")
+    list(APPEND EXTERNAL_LIBRARIES ${P4EST_LIB})
+    list(APPEND EXTERNAL_LIBRARIES ${POPPLER_CPP_LIB})
+    list(APPEND EXTERNAL_INCLUDES ${P4EST_INC})
+  endif()
+endif()
+
 if(HAVE_MESH OR HAVE_SOLVER)
   if(ENABLE_METIS)
     find_library(METIS_LIB metis PATH_SUFFIXES lib)

Common/GmshConfig.h.in

@@ -53,6 +53,7 @@
 #cmakedefine HAVE_OPENGL
 #cmakedefine HAVE_OPTHOM
 #cmakedefine HAVE_OSMESA
+#cmakedefine HAVE_P4EST
 #cmakedefine HAVE_PARSER
 #cmakedefine HAVE_PETSC
 #cmakedefine HAVE_PETSC4PY

Fltk/fieldWindow.cpp

@@ -71,6 +71,7 @@ static void field_put_on_view_cb(Fl_Widget *w, void *data)
 {
   Fl_Menu_Button *mb = ((Fl_Menu_Button *)w);
   Field *field = (Field *)FlGui::instance()->fields->editor_group->user_data();
+  field->update();
   if(mb->value() == 0)
     field->putOnNewView();
   else if(mb->value() - 1 < (int)PView::list.size())

Mesh/automaticMeshSizeField.cpp

 #include "automaticMeshSizeField.h"
+#include "GModel.h"
+#include "GRegion.h"
+#include "MVertex.h"
 
-
+#ifdef HAVE_HXT
+extern "C" {
+#include "hxt_edge.h"
+#include "hxt_curvature.h"
+#include "hxt_bbox.h"
+}
+#endif  
 
 double automaticMeshSizeField:: operator()(double X, double Y, double Z, GEntity *ge) {  
-  return _hbulk;
+#ifdef HAVE_HXT
+  //  return .2;
+  double val = 1.e22;
+  HXTStatus s = hxtOctreeSearchOne(forest, X, Y, Z, &val);
+  if (s == HXT_STATUS_OK){
+    return val;
+  }
+  else Msg::Error ("Cannot find point %g %g %g in the octree",X,Y,Z);
+  return val;
+#else
+  Msg::Error ("Gmsh has to be compiled with HXT for using automaticMeshSizeField");
+#endif  
 }
 
-void automaticMeshSizeField:: update(){
+automaticMeshSizeField::~automaticMeshSizeField(){
+#ifdef HAVE_P4EST
+  if (forest)hxtForestDelete(&forest);  
+  if (forestOptions)hxtForestOptionsDelete(&forestOptions);
+#endif
+}  
+
+
+#ifdef HAVE_HXT
+
+HXTStatus Gmsh2Hxt(std::vector<GRegion *> &regions, HXTMesh *m,
+		   std::map<MVertex *, int> &v2c,
+		   std::vector<MVertex *> &c2v);
+
+
+HXTStatus automaticMeshSizeField:: updateHXT(){
+#ifdef HAVE_P4EST
+
+  
+  //  printf("%d points per circle\n",_nPointsPerCircle);
+  
+  if (!update_needed)return HXT_STATUS_OK;
+  
+  if (forest)hxtForestDelete(&forest);  
+  if (forestOptions)hxtForestOptionsDelete(&forestOptions);
+
+  update_needed = false;
+  // --------------------------------------------------------
+  // Soit on charge un fichier avec l'octree (ma préférence)
+  // Soit on calcule "en live" l'octree ici
+
+  
+  // verify that only one region exists in the model
+  if (GModel::current()->getNumRegions() != 1){
+    Msg::Error ("automaticMeshSizeField can only work for now on models with one region");
+    return HXT_STATUS_ERROR;
+  }
+  std::vector<GRegion*> regions;
+  regions.push_back(*(GModel::current()->firstRegion()));
+  
+  // create HXT mesh structure
+  HXTMesh *mesh;
+  HXTContext *context;
+  HXT_CHECK(hxtContextCreate(&context));
+  HXT_CHECK(hxtMeshCreate(context, &mesh));
+  std::map<MVertex *, int> v2c;
+  std::vector<MVertex *> c2v;  
+  Gmsh2Hxt(regions, mesh, v2c, c2v);
+
+  // Compute curvature
+  double           *curvatureCrossfield;
+  double           *nodalCurvature;
+  HXTEdges         *edges;
+  HXT_CHECK(hxtEdgesCreate(mesh,&edges));
+  HXT_CHECK(hxtCurvatureRusinkiewicz(mesh,&nodalCurvature,&curvatureCrossfield,edges,1));
+  
+  // compute rtree
+  RTree<int , double,3> triRTree;
+  HXTBbox          bbox_triangle;
+
+  for(uint64_t i = 0; i < mesh->triangles.num; ++i){
+    hxtBboxInit(&bbox_triangle);
+    // Noeuds
+    for(int j = 0; j < 3; ++j){
+      double coord[3];
+      int node = mesh->triangles.node[(size_t) 3*i+j];
+      // Coordonnees
+      for(int k = 0; k < 3; ++k){ coord[k] = mesh->vertices.coord[(size_t) 4*node+k]; }
+      hxtBboxAddOne(&bbox_triangle, coord);
+    }
+    int K = (int)i;
+    triRTree.Insert(bbox_triangle.min,bbox_triangle.max,K);
+  }
+  
+  // compute bbox of the mesh
+  double            bbox_vertices[6];
+  HXTBbox           bbox_mesh;
+  hxtBboxInit(&bbox_mesh);
+  // Ajout de tous les points du maillages à la bounding box
+  HXT_CHECK(hxtBboxAdd(&bbox_mesh, mesh->vertices.coord, mesh->vertices.num));
+  for(int i = 0; i < 3; ++i){
+    bbox_vertices[i  ] = 1.3*bbox_mesh.min[i];
+    bbox_vertices[i+3] = 1.3*bbox_mesh.max[i];
+  }
+  
+  // 1 --> OPTIONS ------------------------------------------
+  HXT_CHECK(hxtForestOptionsCreate(&forestOptions));
+  forestOptions->nodePerTwoPi = _nPointsPerCircle;  
+  forestOptions->bbox = bbox_vertices;
+  forestOptions->nodalCurvature = nodalCurvature;
+  forestOptions->triRTree = &triRTree;
+  forestOptions->mesh = mesh;
+  forestOptions->sizeFunction = NULL;
+  
+  // --------------------------------------------------------
+  HXT_CHECK(hxtForestCreate(0, NULL, &forest, NULL, forestOptions));
+  //  return HXT_STATUS_OK;
+  // ---------------------------------------------------------
+  // aller --> niveau 3 par exemple, serait une option ...
+  HXT_CHECK(hxtOctreeRefineToLevel(forest, 3/*_nRefine*/));
+  // idéalement --> convergence
+  int ITER = 0;
+  if (forestOptions->nodePerTwoPi > 0){
+    HXT_CHECK(hxtOctreeCurvatureRefine(forest, 3/*_nRefine*/));
+    // ensuite lissage du gradient 
+    while (ITER++ < 4*_nRefine) {
+      HXT_CHECK(hxtOctreeComputeLaplacian(forest));    
+      HXT_CHECK(hxtOctreeSetMaxGradient(forest));
+    }
+  }
+  forestOptions->nodePerGap = _nPointsPerGap;
+  printf("coucou\n");
+  HXT_CHECK(hxtOctreeSurfacesProches(forest));   
+  ITER = 0;
+  //  while (ITER++ < 4*_nRefine) {
+  //    HXT_CHECK(hxtOctreeComputeLaplacian(forest));    
+  //    HXT_CHECK(hxtOctreeSetMaxGradient(forest));
+  //  }
+
+  HXT_CHECK(hxtEdgesDelete(&edges)       );
+  HXT_CHECK(hxtFree(&curvatureCrossfield));
+  HXT_CHECK(hxtFree(&nodalCurvature)              );
+  HXT_CHECK(hxtMeshDelete(&mesh)                  );
+  HXT_CHECK(hxtContextDelete(&context)            );
+
+
+#endif //ifdef HAVE_P4EST
 }
+#endif
+
+void automaticMeshSizeField:: update(){
+#ifdef HAVE_HXT
+  HXTStatus s = updateHXT();
+  if (s != HXT_STATUS_OK)Msg::Error ("Something went wrong when computing the octree");
+#else
+  Msg::Error ("Gmsh has to be compiled with HXT for using automaticMeshSizeField");
+#endif    
+    
+};

Mesh/automaticMeshSizeField.h

 #ifndef _AUTOMATIC_MESH_SIZE_FIELD_H_
 #define _AUTOMATIC_MESH_SIZE_FIELD_H_
 
+#include "GmshConfig.h"
+
+#ifdef HAVE_HXT
+#include "hxt_octree.h"
+#endif
+
 #include "Field.h"
 
 class automaticMeshSizeField : public Field {
-
+  
+#ifdef HAVE_HXT
+  struct HXTForest *forest;
+  struct HXTForestOptions *forestOptions;
+  HXTStatus updateHXT();
+#endif
+  
   int _nPointsPerCircle;
   int _nPointsPerGap;
   double _hmin, _hmax;
   double _hbulk;
   double _gradientMax;
+  int _nRefine;
   char fileName[256];
-
-public:
-  automaticMeshSizeField()
-  {
-    _nPointsPerCircle = 15;
+  
+ public:
+  ~automaticMeshSizeField();
+  automaticMeshSizeField() : forest(NULL), forestOptions(NULL){
+    _nPointsPerCircle = 55 ;
     _nPointsPerGap = 5;
     _hmin = 1.e-8;// update needed
     _hmax = 1.e+8;// update needed    
     _hbulk = 0.1; // update needed
-    _gradientMax =1.4; 
-
-    
+    _gradientMax =1.4;
+    _nRefine = 5;
     
     options["nPointsPerCircle"] = new FieldOptionInt(_nPointsPerCircle,
-						     "Number of points per circle (adapt to curvature of surfaces)");
+						     "Number of points per circle (adapt to curvature of surfaces)",&update_needed);
     
     options["nPointsPerGap"] = new FieldOptionInt(_nPointsPerGap,
-						  "Number of points in thin layers");
+						  "Number of points in thin layers",&update_needed);
     
     options["hBulk"] = new FieldOptionDouble(_hbulk,
 					     "Size everywhere no size is prescribed", &update_needed);
     
     options["gradientMax"] = new FieldOptionDouble(_gradientMax,
-						   "Maximun gradient of the size field");
+						   "Maximun gradient of the size field",&update_needed);
     
+    options["NRefine"] = new FieldOptionInt(_nRefine,
+					    "Initial refinement level for the octree",&update_needed);
   }
   const char *getName() { return "AutomaticMeshSizeField"; }
   

Mesh/meshGRegionHxt.cpp

@@ -234,9 +234,9 @@ static HXTStatus Hxt2Gmsh(std::vector<GRegion *> &regions, HXTMesh *m,
   return HXT_STATUS_OK;
 }
 
-static HXTStatus Gmsh2Hxt(std::vector<GRegion *> &regions, HXTMesh *m,
-                          std::map<MVertex *, int> &v2c,
-                          std::vector<MVertex *> &c2v)
+HXTStatus Gmsh2Hxt(std::vector<GRegion *> &regions, HXTMesh *m,
+		   std::map<MVertex *, int> &v2c,
+		   std::vector<MVertex *> &c2v)
 {
   std::set<MVertex *> all;
   std::vector<GFace *> faces;

Mesh/meshGRegionHxt.h

@@ -12,4 +12,5 @@ class GRegion;
 
 int meshGRegionHxt(std::vector<GRegion *> &regions);
 
+
 #endif

contrib/hxt/CMakeLists.txt

@@ -39,6 +39,7 @@ set(SRC ${SRC}
   hxt_tetUtils.c
   hxt_tet_aspect_ratio.c
   hxt_vertices.c
+  hxt_octree.cpp
   predicates.c
 )
 endif()

contrib/hxt/hxt_octree.h

+#ifndef HXT_OCTREE_H
+#define HXT_OCTREE_H
+
+// STD INCLUDES
+#include <queue>
+#include <vector>
+
+// GMSH INCLUDES ...
+
+#include "GmshConfig.h"
+#include "SPoint3.h"
+#include "SVector3.h"
+#include "rtree.h"
+
+// P4EST INCLUDES
+
+#ifdef HAVE_P4EST
+#include <p4est_to_p8est.h>
+#include <p8est_bits.h>
+#include <p8est_ghost.h>
+#include <p8est_nodes.h>
+#include <p8est_vtk.h>
+#include <p8est_iterate.h>
+#include <p8est_extended.h>
+#include <p8est_search.h>
+#endif 
+
+// HXT INCLUDES
+
+extern "C" {
+#include "hxt_api.h"
+#include "hxt_mesh.h"
+#include "hxt_bbox.h"
+}
+
+typedef struct HXTForestOptions{
+  int                   nodePerTwoPi;
+  int                   nodePerGap;
+  double               *bbox;
+  double               *nodalCurvature;
+  double              (*sizeFunction)(double, double, double) ;
+  RTree<int,double,3>  *triRTree;
+  HXTMesh              *mesh;
+} HXTForestOptions;
+
+typedef struct HXTForest{
+#ifdef HAVE_P4EST
+  p4est_t          *p4est;
+#endif 
+  HXTForestOptions  *forestOptions;
+} HXTForest;
+
+// Donnees disponibles sur chaque quadrant
+typedef struct size_data{
+  double size;
+#ifdef HAVE_P4EST
+  double ds[P4EST_DIM];
+#endif 
+  double d2s; // Laplacien
+
+  double h;
+  // Les tailles pour les differences finies
+  double h_xL, h_xR;
+  double h_yD, h_yU;
+  
+  double h_xavg;
+  double h_yavg;
+  
+  double h_zB, h_zT;
+  double h_zavg;
+
+  double hMin;
+  
+  // Pour raffiner sur la courbure : pas possible de donner un user_pointer
+  //  à p4est_refine ?
+  int refineFlag; 
+  
+} size_data_t;
+
+// Un point pour la recherche dans l'arbre
+typedef struct size_point{
+  double x;
+  double y;
+  double z;
+  double size;
+} size_point_t;
+
+typedef struct size_fun{
+  double (*myFun)(double, double, double);
+} size_fun_t;
+
+// API ---------------------------------------------------------------------------------------------
+
+HXTStatus hxtOctreeRefineToLevel(HXTForest *forest, int lvl);
+HXTStatus hxtOctreeComputeLaplacian(HXTForest *forest);
+HXTStatus hxtOctreeLaplacianRefine(HXTForest *forest, int nRefine);
+HXTStatus hxtOctreeSetMaxGradient(HXTForest *forest);
+HXTStatus hxtOctreeSmoothGradient(HXTForest *forest, int nMax);
+HXTStatus hxtForestOptionsCreate(HXTForestOptions **forestOptions);
+HXTStatus hxtForestOptionsDelete(HXTForestOptions **forestOptions);
+HXTStatus hxtForestCreate(int argc, char **argv, HXTForest **forest, const char* filename, HXTForestOptions *forestOptions);
+HXTStatus hxtForestDelete(HXTForest **forest);
+HXTStatus hxtOctreeRTreeIntersection(HXTForest *forest);
+HXTStatus hxtOctreeCurvatureRefine(HXTForest *forest, int nMax);
+HXTStatus hxtOctreeSearchOne(HXTForest *forest, double x, double y, double z, double *size);
+HXTStatus hxtOctreeSurfacesProches(HXTForest *forest);
+HXTStatus hxtOctreeElementEstimation(HXTForest *forest, double *elemEstimate);
+
+#endif