Select Git revision
automaticMeshSizeField.cpp
automaticMeshSizeField.cpp 5.04 KiB
#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)
{
double val = 1.e22;
#if defined(HAVE_HXT) && defined(HAVE_P4EST)
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);
#else
Msg::Error ("Gmsh has to be compiled with HXT and P4EST for using automaticMeshSizeField");
#endif
return val;
}
automaticMeshSizeField::~automaticMeshSizeField(){
#if defined(HAVE_HXT) && defined(HAVE_P4EST)
if (forest)hxtForestDelete(&forest);
if (forestOptions)hxtForestOptionsDelete(&forestOptions);
#endif
}
#if defined(HAVE_HXT) && defined(HAVE_P4EST)
HXTStatus Gmsh2Hxt(std::vector<GRegion *> ®ions, HXTMesh *m,
std::map<MVertex *, int> &v2c,
std::vector<MVertex *> &c2v);
HXTStatus automaticMeshSizeField:: updateHXT(){
// 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
void automaticMeshSizeField:: update(){
#if defined(HAVE_HXT) && defined(HAVE_P4EST)
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 and P4EST for using automaticMeshSizeField");
#endif
};