From 49bbacfcbfda3571c3e7ef9440df68bb9f64dabc Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Tue, 8 Jul 2008 12:43:26 +0000 Subject: [PATCH] added octree to GModel for fast element lookups --- Common/Octree.cpp | 53 +++++++++--------- Common/Octree.h | 4 +- Common/OctreeInternals.h | 5 +- Geo/GModel.cpp | 84 +++++++++++++++++++++++++++-- Geo/GModel.h | 12 +++++ Mesh/meshGEdge.cpp | 3 +- Mesh/meshGFace.cpp | 3 +- Mesh/meshGRegion.cpp | 3 +- Post/OctreePost.cpp | 114 ++++++++++++--------------------------- Post/OctreePost.h | 1 - Post/PViewDataGModel.h | 2 + 11 files changed, 163 insertions(+), 121 deletions(-) diff --git a/Common/Octree.cpp b/Common/Octree.cpp index 9dfbd4b7e5..231de3a280 100644 --- a/Common/Octree.cpp +++ b/Common/Octree.cpp @@ -1,4 +1,4 @@ -// $Id: Octree.cpp,v 1.5 2008-03-20 11:44:02 geuzaine Exp $ +// $Id: Octree.cpp,v 1.6 2008-07-08 12:43:25 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -24,10 +24,6 @@ #include <list> #include "Octree.h" -using std::list; - -void free_buckets(octantBucket *); - Octree* Octree_Create(int maxElements, double origin[3], double size[3], void (*BB)(void *, double*, double*), void (*Centroid)(void *, double *), @@ -42,11 +38,34 @@ Octree* Octree_Create(int maxElements, double origin[3], double size[3], return myOctree; } +void free_buckets(octantBucket * bucket) +{ + int i, numBuck = 8; + ELink ptr1, ptr2; + + if(bucket->next == NULL) { + ptr1 = bucket->lhead; + while(ptr1 != NULL) { + ptr2 = ptr1; + ptr1 = ptr1->next; + delete ptr2; + } + bucket->listBB.clear(); + return; + } + + for(i = numBuck-1; i >= 0; i--) + free_buckets((bucket->next)+i); + delete []bucket->next; + return; +} + void Octree_Delete(Octree *myOctree) { delete myOctree->info; free_buckets(myOctree->root); delete myOctree->root; + delete myOctree; } void Octree_Insert(void * element, Octree *myOctree) @@ -78,29 +97,7 @@ void * Octree_Search(double *pt, Octree *myOctree) myOctree->function_BB, myOctree->function_inElement); } -void free_buckets(octantBucket * bucket) -{ - int i, numBuck = 8; - ELink ptr1, ptr2; - - if(bucket->next == NULL) { - ptr1 = bucket->lhead; - while(ptr1 != NULL) { - ptr2 = ptr1; - ptr1 = ptr1->next; - delete ptr2; - } - bucket->listBB.clear(); - return; - } - - for(i = numBuck-1; i >= 0; i--) - free_buckets((bucket->next)+i); - delete []bucket->next; - return; -} - -void Octree_SearchAll(double * pt, Octree * myOctree, list<void *> * output) +void Octree_SearchAll(double * pt, Octree * myOctree, std::list<void *> * output) { searchAllElements(myOctree->root, pt, myOctree->info, myOctree->function_BB, myOctree->function_inElement, output); diff --git a/Common/Octree.h b/Common/Octree.h index 92a28b93b9..ce3aa396d4 100644 --- a/Common/Octree.h +++ b/Common/Octree.h @@ -23,8 +23,6 @@ #include <list> #include "OctreeInternals.h" -using std::list; - Octree* Octree_Create(int maxElements, // max. num of elts allowed in an octant double *origin, // smallest x,y, z of model's bounding box double *size, // size in x, y, z of model bounding box @@ -36,6 +34,6 @@ void Octree_Delete(Octree *); void Octree_Insert(void *, Octree *); void Octree_Arrange(Octree *); void * Octree_Search(double *, Octree *); -void Octree_SearchAll(double *, Octree *, list<void *> *); +void Octree_SearchAll(double *, Octree *, std::list<void *> *); #endif diff --git a/Common/OctreeInternals.h b/Common/OctreeInternals.h index cf90792f83..7037cf6c4c 100644 --- a/Common/OctreeInternals.h +++ b/Common/OctreeInternals.h @@ -62,14 +62,15 @@ struct global { }; typedef struct global globalInfo; -typedef struct +class Octree { + public: globalInfo *info; octantBucket *root; BBFunction function_BB; InEleFunction function_inElement; CentroidFunction function_centroid; -}Octree; +}; void refineOctants(octantBucket *buckets, globalInfo *globalPara); diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp index 6435db5733..33995d971e 100644 --- a/Geo/GModel.cpp +++ b/Geo/GModel.cpp @@ -1,4 +1,4 @@ -// $Id: GModel.cpp,v 1.94 2008-07-03 18:15:29 geuzaine Exp $ +// $Id: GModel.cpp,v 1.95 2008-07-08 12:43:25 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -32,6 +32,7 @@ #else # include "Message.h" # include "gmshSurface.h" +# include "Octree.h" # include "Field.h" # include "Generator.h" # include "Context.h" @@ -43,8 +44,8 @@ std::vector<GModel*> GModel::list; int GModel::_current = -1; GModel::GModel(std::string name) - : _geo_internals(0), _occ_internals(0), _fm_internals(0), _fields(0), - _currentMeshEntity(0), modelName(name), normals(0) + : _octree(0), _geo_internals(0), _occ_internals(0), _fm_internals(0), + _fields(0), _currentMeshEntity(0), modelName(name), normals(0) { list.push_back(this); @@ -115,8 +116,7 @@ void GModel::destroy() if(normals) delete normals; normals = 0; - _vertexVectorCache.clear(); - _vertexMapCache.clear(); + destroyMeshCaches(); MVertex::resetGlobalNumber(); MElement::resetGlobalNumber(); @@ -127,6 +127,14 @@ void GModel::destroy() #endif } +void GModel::destroyMeshCaches() +{ + _vertexVectorCache.clear(); + _vertexMapCache.clear(); + if(_octree) Octree_Delete(_octree); + _octree = 0; +} + GRegion *GModel::getRegionByTag(int n) const { GEntity tmp((GModel*)this, n); @@ -443,6 +451,72 @@ int GModel::getNumMeshElements() return n; } +static void MElementBB(void *a, double *min, double *max) +{ + MElement *e = (MElement*)a; + MVertex *v = e->getVertex(0); + min[0] = max[0] = v->x(); + min[1] = max[1] = v->y(); + min[2] = max[2] = v->z(); + for(int i = 1; i < e->getNumVertices(); i++){ + v = e->getVertex(i); + min[0] = std::min(min[0], v->x()); max[0] = std::max(max[0], v->x()); + min[1] = std::min(min[1], v->y()); max[1] = std::max(max[1], v->y()); + min[2] = std::min(min[2], v->z()); max[2] = std::max(max[2], v->z()); + } +} + +static void MElementCentroid(void *a, double *x) +{ + MElement *e = (MElement*)a; + MVertex *v = e->getVertex(0); + int n = e->getNumVertices(); + x[0] = v->x(); x[1] = v->y(); x[2] = v->z(); + for(int i = 1; i < n; i++) { + v = e->getVertex(i); + x[0] += v->x(); x[1] += v->y(); x[2] += v->z(); + } + double oc = 1. / (double)n; + x[0] *= oc; + x[1] *= oc; + x[2] *= oc; +} + +static int MElementInEle(void *a, double *x) +{ + MElement *e = (MElement*)a; + double uvw[3]; + e->xyz2uvw(x, uvw); + return e->isInside(uvw[0], uvw[1], uvw[2]) ? 1 : 0; +} + +MElement *GModel::getMeshElementByCoord(SPoint3 &p) +{ +#if !defined(HAVE_GMSH_EMBEDDED) + if(!_octree){ + Msg::Debug("Rebuilding mesh element octree"); + SBoundingBox3d bb = bounds(); + double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()}; + double size[3] = {bb.max().x() - bb.min().x(), + bb.max().y() - bb.min().y(), + bb.max().z() - bb.min().z()}; + const int maxElePerBucket = 100; // memory vs. speed trade-off + _octree = Octree_Create(maxElePerBucket, min, size, + MElementBB, MElementCentroid, MElementInEle); + std::vector<GEntity*> ent = getEntities(); + for(unsigned int i = 0; i < ent.size(); i++) + for(unsigned int j = 0; j < ent[i]->getNumMeshElements(); j++) + Octree_Insert(ent[i]->getMeshElement(j), _octree); + Octree_Arrange(_octree); + } + double P[3] = {p.x(), p.y(), p.z()}; + return (MElement*)Octree_Search(P, _octree); +#else + Msg::Error("Embedded Gmsh cannot perform octree-based element searches"); + return 0; +#endif +} + template <class T> static void insertMeshVertices(std::vector<MVertex*> &vertices, T &container) { diff --git a/Geo/GModel.h b/Geo/GModel.h index d7d76f35a6..0d3bb6bf74 100644 --- a/Geo/GModel.h +++ b/Geo/GModel.h @@ -27,8 +27,10 @@ #include "GEdge.h" #include "GFace.h" #include "GRegion.h" +#include "SPoint3.h" #include "SBoundingBox3d.h" +class Octree; class FM_Internals; class GEO_Internals; class OCC_Internals; @@ -44,6 +46,9 @@ class GModel std::vector<MVertex*> _vertexVectorCache; std::map<int, MVertex*> _vertexMapCache; + // an octree for fast mesh element lookup + Octree *_octree; + // geo model internal data GEO_Internals *_geo_internals; void _createGEOInternals(); @@ -103,6 +108,10 @@ class GModel // delete everything in a GModel void destroy(); + // delete all the mesh-related caches (this must be called when the + // mesh is changed) + void destroyMeshCaches(); + // access internal CAD representations GEO_Internals *getGEOInternals(){ return _geo_internals; } OCC_Internals *getOCCInternals(){ return _occ_internals; } @@ -196,6 +205,9 @@ class GModel // Returns the total number of elements in the mesh int getNumMeshElements(); + // Access a mesh element by coordinates + MElement *getMeshElementByCoord(SPoint3 &p); + // Returns the total number of vertices in the mesh int getNumMeshVertices(); diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp index 15729f03ca..41ac4d298f 100644 --- a/Mesh/meshGEdge.cpp +++ b/Mesh/meshGEdge.cpp @@ -1,4 +1,4 @@ -// $Id: meshGEdge.cpp,v 1.67 2008-07-03 18:15:29 geuzaine Exp $ +// $Id: meshGEdge.cpp,v 1.68 2008-07-08 12:43:25 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -275,6 +275,7 @@ void deMeshGEdge::operator() (GEdge *ge) delete ge->lines[i]; ge->lines.clear(); ge->deleteVertexArrays(); + ge->model()->destroyMeshCaches(); } void meshGEdge::operator() (GEdge *ge) diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index 0d50915425..4976023ecd 100644 --- a/Mesh/meshGFace.cpp +++ b/Mesh/meshGFace.cpp @@ -1,4 +1,4 @@ -// $Id: meshGFace.cpp,v 1.140 2008-07-03 18:15:29 geuzaine Exp $ +// $Id: meshGFace.cpp,v 1.141 2008-07-08 12:43:25 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -1292,6 +1292,7 @@ void deMeshGFace::operator() (GFace *gf) for (unsigned int i=0;i<gf->quadrangles.size();i++) delete gf->quadrangles[i]; gf->quadrangles.clear(); gf->deleteVertexArrays(); + gf->model()->destroyMeshCaches(); gf->meshStatistics.status = GFace::PENDING; gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0; diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp index 81511767d3..db73d425fe 100644 --- a/Mesh/meshGRegion.cpp +++ b/Mesh/meshGRegion.cpp @@ -1,4 +1,4 @@ -// $Id: meshGRegion.cpp,v 1.56 2008-07-03 18:15:29 geuzaine Exp $ +// $Id: meshGRegion.cpp,v 1.57 2008-07-08 12:43:26 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -428,6 +428,7 @@ void deMeshGRegion::operator() (GRegion *gr) delete gr->pyramids[i]; gr->pyramids.clear(); gr->deleteVertexArrays(); + gr->model()->destroyMeshCaches(); } int intersect_line_triangle(double X[3], double Y[3], double Z[3] , diff --git a/Post/OctreePost.cpp b/Post/OctreePost.cpp index 8e633cacbd..ca2ca837d3 100644 --- a/Post/OctreePost.cpp +++ b/Post/OctreePost.cpp @@ -1,4 +1,4 @@ -// $Id: OctreePost.cpp,v 1.15 2008-06-07 17:20:57 geuzaine Exp $ +// $Id: OctreePost.cpp,v 1.16 2008-07-08 12:43:26 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -219,47 +219,6 @@ static void addListOfStuff(Octree *o, List_T *l, int nbelm) } } -// helper routines for GModel-based views - -void MElementBB(void *a, double *min, double *max) -{ - MElement *e = (MElement*)a; - MVertex *v = e->getVertex(0); - min[0] = max[0] = v->x(); - min[1] = max[1] = v->y(); - min[2] = max[2] = v->z(); - for(int i = 1; i < e->getNumVertices(); i++){ - v = e->getVertex(i); - min[0] = std::min(min[0], v->x()); max[0] = std::max(max[0], v->x()); - min[1] = std::min(min[1], v->y()); max[1] = std::max(max[1], v->y()); - min[2] = std::min(min[2], v->z()); max[2] = std::max(max[2], v->z()); - } -} - -void MElementCentroid(void *a, double *x) -{ - MElement *e = (MElement*)a; - MVertex *v = e->getVertex(0); - int n = e->getNumVertices(); - x[0] = v->x(); x[1] = v->y(); x[2] = v->z(); - for(int i = 1; i < n; i++) { - v = e->getVertex(i); - x[0] += v->x(); x[1] += v->y(); x[2] += v->z(); - } - double oc = 1. / (double)n; - x[0] *= oc; - x[1] *= oc; - x[2] *= oc; -} - -int MElementInEle(void *a, double *x) -{ - MElement *e = (MElement*)a; - double uvw[3]; - e->xyz2uvw(x, uvw); - return e->isInside(uvw[0], uvw[1], uvw[2]) ? 1 : 0; -} - // OctreePost implementation OctreePost::~OctreePost() @@ -270,26 +229,27 @@ OctreePost::~OctreePost() if(_SH) Octree_Delete(_SH); if(_VH) Octree_Delete(_VH); if(_TH) Octree_Delete(_TH); if(_SI) Octree_Delete(_SI); if(_VI) Octree_Delete(_VI); if(_TI) Octree_Delete(_TI); if(_SY) Octree_Delete(_SY); if(_VY) Octree_Delete(_VY); if(_TY) Octree_Delete(_TY); - if(_GModel) Octree_Delete(_GModel); } OctreePost::OctreePost(PView *v) : _ST(0), _VT(0), _TT(0), _SQ(0), _VQ(0), _TQ(0), _SS(0), _VS(0), _TS(0), _SH(0), _VH(0), _TH(0), _SI(0), _VI(0), _TI(0), _SY(0), _VY(0), _TY(0), - _GModel(0), _theView(v), _theViewDataList(0), _theViewDataGModel(0) + _theView(v), _theViewDataList(0), _theViewDataGModel(0) { - SBoundingBox3d bb = v->getData()->getBoundingBox(); - - double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()}; - - double size[3] = {bb.max().x() - bb.min().x(), - bb.max().y() - bb.min().y(), - bb.max().z() - bb.min().z()}; + _theViewDataGModel = dynamic_cast<PViewDataGModel*>(_theView->getData()); - const int maxElePerBucket = 100; // memory vs. speed trade-off + if(_theViewDataGModel) return; // the octree is available in the model _theViewDataList = dynamic_cast<PViewDataList*>(_theView->getData()); + if(_theViewDataList){ + SBoundingBox3d bb = v->getData()->getBoundingBox(); + double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()}; + double size[3] = {bb.max().x() - bb.min().x(), + bb.max().y() - bb.min().y(), + bb.max().z() - bb.min().z()}; + const int maxElePerBucket = 100; // memory vs. speed trade-off + PViewDataList *l = _theViewDataList; _SL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle); addListOfStuff(_SL, l->SL, 6 + 2 * l->getNumTimeSteps()); @@ -300,7 +260,7 @@ OctreePost::OctreePost(PView *v) _TL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle); addListOfStuff(_TL, l->TL, 6 + 18 * l->getNumTimeSteps()); Octree_Arrange(_TL); - + _ST = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle); addListOfStuff(_ST, l->ST, 9 + 3 * l->getNumTimeSteps()); Octree_Arrange(_ST); @@ -310,7 +270,7 @@ OctreePost::OctreePost(PView *v) _TT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle); addListOfStuff(_TT, l->TT, 9 + 27 * l->getNumTimeSteps()); Octree_Arrange(_TT); - + _SQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle); addListOfStuff(_SQ, l->SQ, 12 + 4 * l->getNumTimeSteps()); Octree_Arrange(_SQ); @@ -320,7 +280,7 @@ OctreePost::OctreePost(PView *v) _TQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle); addListOfStuff(_TQ, l->TQ, 12 + 36 * l->getNumTimeSteps()); Octree_Arrange(_TQ); - + _SS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle); addListOfStuff(_SS, l->SS, 12 + 4 * l->getNumTimeSteps()); Octree_Arrange(_SS); @@ -330,7 +290,7 @@ OctreePost::OctreePost(PView *v) _TS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle); addListOfStuff(_TS, l->TS, 12 + 36 * l->getNumTimeSteps()); Octree_Arrange(_TS); - + _SH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle); addListOfStuff(_SH, l->SH, 24 + 8 * l->getNumTimeSteps()); Octree_Arrange(_SH); @@ -350,7 +310,7 @@ OctreePost::OctreePost(PView *v) _TI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle); addListOfStuff(_TI, l->TI, 18 + 54 * l->getNumTimeSteps()); Octree_Arrange(_TI); - + _SY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle); addListOfStuff(_SY, l->SY, 15 + 5 * l->getNumTimeSteps()); Octree_Arrange(_SY); @@ -361,22 +321,6 @@ OctreePost::OctreePost(PView *v) addListOfStuff(_TY, l->TY, 15 + 45 * l->getNumTimeSteps()); Octree_Arrange(_TY); } - else{ - _theViewDataGModel = dynamic_cast<PViewDataGModel*>(_theView->getData()); - if(_theViewDataGModel){ - // we should think about storing the octree in the model, so - // that we can reuse it multiple times - _GModel = Octree_Create(maxElePerBucket, min, size, - MElementBB, MElementCentroid, MElementInEle); - for(int i = 0; i < _theViewDataGModel->getNumEntities(0); i++){ - GEntity *ge = _theViewDataGModel->getEntity(0, i); - for(unsigned int j = 0; j < ge->getNumMeshElements(); j++) - Octree_Insert(ge->getMeshElement(j), _GModel); - } - Octree_Arrange(_GModel); - } - } - } bool OctreePost::_getValue(void *in, int dim, int nbNod, int nbComp, @@ -410,8 +354,8 @@ bool OctreePost::_getValue(void *in, int dim, int nbNod, int nbComp, return true; } -bool OctreePost::_getValue(void *in, int nbComp, double P[3], int timestep, double *values, - double *elementSize) +bool OctreePost::_getValue(void *in, int nbComp, double P[3], int timestep, + double *values, double *elementSize) { if(!in) return false; @@ -475,7 +419,11 @@ bool OctreePost::searchScalar(double x, double y, double z, double *values, if(_getValue(Octree_Search(P, _SL), 1, 2, 1, P, step, values, size)) return true; } else if(_theViewDataGModel){ - if(_getValue(Octree_Search(P, _GModel), 1, P, step, values, size)) return true; + GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step); + if(m){ + SPoint3 pt(P); + if(_getValue(m->getMeshElementByCoord(pt), 1, P, step, values, size)) return true; + } } return false; @@ -503,9 +451,13 @@ bool OctreePost::searchVector(double x, double y, double z, double *values, if(_getValue(Octree_Search(P, _VL), 1, 2, 3, P, step, values, size)) return true; } else if(_theViewDataGModel){ - if(_getValue(Octree_Search(P, _GModel), 3, P, step, values, size)) return true; + GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step); + if(m){ + SPoint3 pt(P); + if(_getValue(m->getMeshElementByCoord(pt), 3, P, step, values, size)) return true; + } } - + return false; } @@ -531,7 +483,11 @@ bool OctreePost::searchTensor(double x, double y, double z, double *values, if(_getValue(Octree_Search(P, _TL), 1, 2, 9, P, step, values, size)) return true; } else if(_theViewDataGModel){ - if(_getValue(Octree_Search(P, _GModel), 9, P, step, values, size)) return true; + GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step); + if(m){ + SPoint3 pt(P); + if(_getValue(m->getMeshElementByCoord(pt), 9, P, step, values, size)) return true; + } } return false; diff --git a/Post/OctreePost.h b/Post/OctreePost.h index 110edff831..69ed7d2016 100644 --- a/Post/OctreePost.h +++ b/Post/OctreePost.h @@ -36,7 +36,6 @@ class OctreePost Octree *_SH, *_VH, *_TH; Octree *_SI, *_VI, *_TI; Octree *_SY, *_VY, *_TY; - Octree *_GModel; PView *_theView; PViewDataList *_theViewDataList; PViewDataGModel *_theViewDataGModel; diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h index 30189e325d..7f651a9dcd 100644 --- a/Post/PViewDataGModel.h +++ b/Post/PViewDataGModel.h @@ -186,6 +186,8 @@ class PViewDataGModel : public PViewData { GEntity *getEntity(int step, int ent); // direct access to value by index bool getValueByIndex(int step, int dataIndex, int node, int comp, double &val); + // get underlying model + GModel* getModel(int step){ return _steps[step]->getModel(); } // Add some data "on the fly" (data is stored by vertex: if a field // has e.g. 3 components, nodalData contains 3 * N entries with N -- GitLab