diff --git a/CMakeLists.txt b/CMakeLists.txt index 9852f327445212e2022931993aea5f2554bf15b5..d23017c98d70f1501841195ce9dfa2a7195936cd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,7 +81,6 @@ opt(PLUGINS "Enable post-processing plugins" ${DEFAULT}) opt(POST "Enable post-processing module (required by GUI)" ${DEFAULT}) opt(POPPLER "Enable Poppler for displaying PDF documents (experimental)" OFF) opt(QT "Enable dummy QT graphical interface proof-of-concept (experimental)" OFF) -opt(RTREE "Enable RTREE (used for quad/hex mesh generation)" ${DEFAULT}) opt(SALOME "Enable Salome routines for CAD healing" ${DEFAULT}) opt(SGEOM "Enable SGEOM interface to OCC (experimental)" OFF) opt(SLEPC "Enable SLEPc eigensolvers (required for conformal compounds)" ${DEFAULT}) @@ -756,11 +755,6 @@ if(HAVE_MESH OR HAVE_SOLVER) endif(HAVE_MESH OR HAVE_SOLVER) if(HAVE_MESH) - if(ENABLE_RTREE) - include_directories(contrib/rtree) - set_config_option(HAVE_RTREE "RTree") - endif(ENABLE_RTREE) - if(ENABLE_VORO3D) add_subdirectory(contrib/voro++) include_directories(contrib/voro++/src) diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in index 98118893fd231f03afefc3a7b4ae2503fd6a88d8..c2da4a4d345284c379b800ad4540d9fcc6b19aff 100644 --- a/Common/GmshConfig.h.in +++ b/Common/GmshConfig.h.in @@ -60,7 +60,6 @@ #cmakedefine HAVE_POST #cmakedefine HAVE_POPPLER #cmakedefine HAVE_QT -#cmakedefine HAVE_RTREE #cmakedefine HAVE_SALOME #cmakedefine HAVE_SGEOM #cmakedefine HAVE_SLEPC diff --git a/contrib/rtree/rtree.h b/Common/rtree.h similarity index 97% rename from contrib/rtree/rtree.h rename to Common/rtree.h index 31f21e579a62d738e6ac506dc7edf7c03f1e8e2f..c4ea149046a5b9a12dbefd0b94e787bd757177f1 100644 --- a/contrib/rtree/rtree.h +++ b/Common/rtree.h @@ -1,3 +1,49 @@ +/* + +TITLE + + R-TREES: A DYNAMIC INDEX STRUCTURE FOR SPATIAL SEARCHING + +DESCRIPTION + + A C++ templated version of the RTree algorithm. + For more information please read the comments in RTree.h + +AUTHORS + + * 1983 Original algorithm and test code by Antonin Guttman and Michael Stonebraker, UC Berkely + * 1994 ANCI C ported from original test code by Melinda Green - melinda@superliminal.com + * 1995 Sphere volume fix for degeneracy problem submitted by Paul Brook + * 2004 Templated C++ port by Greg Douglas + +LICENSE: + + Entirely free for all uses. Enjoy! + +FILES + * RTree.h The C++ templated RTree implementation. Well commented. + * Test.cpp A simple test program, ported from the original C version. + * MemoryTest.cpp A more rigourous test to validate memory use. + * README.TXT This file. + +TO BUILD + + To build a test, compile only one of the test files with RTree.h. + Both test files contain a main() function. + +RECENT CHANGE LOG + +05 Jan 2010 +o Fixed Iterator GetFirst() - Previous fix was not incomplete + +03 Dec 2009 +o Added Iteartor GetBounds() +o Added Iterator usage to simple test +o Fixed Iterator GetFirst() - Thanks Mathew Riek +o Minor updates for MSVC 2005/08 compilers + + */ + #ifndef RTREE_H #define RTREE_H #include <algorithm> diff --git a/Mesh/filterElements.cpp b/Mesh/filterElements.cpp index fa25fc70bc2db37265d48e1786b6a330295af99f..41fd2d0e0db4d04cb13569acf124ae619a3a61df 100644 --- a/Mesh/filterElements.cpp +++ b/Mesh/filterElements.cpp @@ -15,9 +15,8 @@ #include "MQuadrangle.h" #include "MPrism.h" #include "MHexahedron.h" - -#if defined(HAVE_RTREE) #include "rtree.h" + void MElementBB(void *a, double *min, double *max); int MElementInEle(void *a, double *x); @@ -221,29 +220,3 @@ void filterOverlappingElements (std::vector<MElement*> &els, els = newEls; } -#else - -void filterOverlappingElements (std::vector<MTriangle*> &blTris, - std::vector<MQuadrangle*>&blQuads, - std::map<MElement*,std::vector<MElement*> > &_elemColumns, - std::map<MElement*,MElement*> &_toFirst) -{ - Msg::Error("Gmsh needs to be compiled with RTREE support for bonudary layers"); -} - -void filterOverlappingElements (std::vector<MPrism*> &blPrisms, - std::vector<MHexahedron*>&blHexes, - std::map<MElement*,std::vector<MElement*> > &_elemColumns, - std::map<MElement*,MElement*> &_toFirst) -{ - Msg::Error("Gmsh needs to be compiled with RTREE support for bonudary layers"); -} - -void filterOverlappingElements (std::vector<MElement*> &els, - std::map<MElement*,std::vector<MElement*> > &_elemColumns, - std::map<MElement*,MElement*> &_toFirst ) -{ - Msg::Error("Gmsh needs to be compiled with RTREE support for bonudary layers"); -} - -#endif diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp index 1b887c40359520f0d50def2d86a1f9355b546ac0..86a734c8c5b12f9a4030dc11dd6dd8161be79954 100644 --- a/Mesh/meshGFaceExtruded.cpp +++ b/Mesh/meshGFaceExtruded.cpp @@ -55,7 +55,7 @@ static void createQuaTri(std::vector<MVertex*> &v, GFace *to, else addQuadrangle(v[0], v[1], v[3], v[2], to); } - + } else if(!constrainedEdges){ addTriangle(v[0], v[1], v[3], to); @@ -179,23 +179,24 @@ static void copyMesh(GFace *from, GFace *to, int quadToTri_valid = IsValidQuadToTriTop(to, &quadToTri, &detectQuadToTriTop); bool is_toroidal = quadToTri_valid >= 2 ? true : false; bool is_noaddverts = quadToTri_valid == 3 ? true : false; - if( detectQuadToTriTop && !quadToTri_valid && !is_toroidal ){ + if(detectQuadToTriTop && !quadToTri_valid && !is_toroidal){ Msg::Error("In MeshGFaceExtrudedSurface::copyMesh(), Mesh of QuadToTri top " "surface %d likely has errors.", to->tag()); } - - // if this is toroidal No New Vertices QuadToTri, then replace the root dependency face's boundary - // quads with triangles for better meshing. - if( is_toroidal && is_noaddverts ){ - GFace *root = findRootSourceFaceForFace( from ); - if( root == from ){ - ReplaceBndQuadsInFace( root ); - Msg::Warning("To facilitate QuadToTri interface on surface %d, source surface %d was re-meshed " - "with all triangles on boundary. To avoid this, use QuadTriAddVerts instead of " - "QuadTriNoNewVerts.", to->tag(), root->tag()); + + // if this is toroidal No New Vertices QuadToTri, then replace the root + // dependency face's boundary quads with triangles for better meshing. + if(is_toroidal && is_noaddverts){ + GFace *root = findRootSourceFaceForFace(from); + if(root == from){ + ReplaceBndQuadsInFace(root); + Msg::Warning("To facilitate QuadToTri interface on surface %d, source " + "surface %d was re-meshed with all triangles on boundary. " + "To avoid this, use QuadTriAddVerts instead of QuadTriNoNewVerts", + to->tag(), root->tag()); } } - + // create triangle elements std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp; for(unsigned int i = 0; i < from->triangles.size(); i++){ @@ -227,7 +228,7 @@ static void copyMesh(GFace *from, GFace *to, Msg::Error("In MeshExtrudedSurface()::copyMesh(), mesh of QuadToTri top " "surface %d failed.", to->tag() ); return; - } + } // create quadrangle elements if NOT QuadToTri and NOT toroidal diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp index d896aa88b7c11bd76f8e9c665b6cf6eeda32a31a..16eaf94fd4347ec7b953de931a32de4f75f748ab 100644 --- a/Mesh/simple3D.cpp +++ b/Mesh/simple3D.cpp @@ -19,10 +19,7 @@ #include "Context.h" #include <iostream> #include <string> - -#if defined(HAVE_RTREE) #include "rtree.h" -#endif #define k1 0.7 //k1*h is the minimal distance between two nodes #define k2 0.5 //k2*h is the minimal distance to the boundary @@ -100,7 +97,8 @@ class Wrapper{ /*********functions*********/ -double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m){ +double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m) +{ double distance; double x1,y1,z1; double x2,y2,z2; @@ -128,7 +126,8 @@ double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m){ return distance; } -bool rtree_callback(Node* neighbour,void* w){ +bool rtree_callback(Node* neighbour,void* w) +{ double h; double distance; Metric m; @@ -154,7 +153,8 @@ bool rtree_callback(Node* neighbour,void* w){ /*********class Metric*********/ -Metric::Metric(){ +Metric::Metric() +{ m11 = 1.0; m21 = 0.0; m31 = 0.0; @@ -168,77 +168,41 @@ Metric::Metric(){ Metric::~Metric(){} -void Metric::set_m11(double new_m11){ - m11 = new_m11; -} +void Metric::set_m11(double new_m11){ m11 = new_m11; } -void Metric::set_m21(double new_m21){ - m21 = new_m21; -} +void Metric::set_m21(double new_m21){ m21 = new_m21; } -void Metric::set_m31(double new_m31){ - m31 = new_m31; -} +void Metric::set_m31(double new_m31){ m31 = new_m31; } -void Metric::set_m12(double new_m12){ - m12 = new_m12; -} +void Metric::set_m12(double new_m12){ m12 = new_m12; } -void Metric::set_m22(double new_m22){ - m22 = new_m22; -} +void Metric::set_m22(double new_m22){ m22 = new_m22; } -void Metric::set_m32(double new_m32){ - m32 = new_m32; -} +void Metric::set_m32(double new_m32){ m32 = new_m32; } -void Metric::set_m13(double new_m13){ - m13 = new_m13; -} +void Metric::set_m13(double new_m13){ m13 = new_m13; } -void Metric::set_m23(double new_m23){ - m23 = new_m23; -} +void Metric::set_m23(double new_m23){ m23 = new_m23; } -void Metric::set_m33(double new_m33){ - m33 = new_m33; -} +void Metric::set_m33(double new_m33){ m33 = new_m33; } -double Metric::get_m11(){ - return m11; -} +double Metric::get_m11(){ return m11; } -double Metric::get_m21(){ - return m21; -} +double Metric::get_m21(){ return m21; } -double Metric::get_m31(){ - return m31; -} +double Metric::get_m31(){ return m31; } -double Metric::get_m12(){ - return m12; -} +double Metric::get_m12(){ return m12; } -double Metric::get_m22(){ - return m22; -} +double Metric::get_m22(){ return m22; } -double Metric::get_m32(){ - return m32; -} +double Metric::get_m32(){ return m32; } -double Metric::get_m13(){ - return m13; -} +double Metric::get_m13(){ return m13; } -double Metric::get_m23(){ - return m23; -} +double Metric::get_m23(){ return m23; } -double Metric::get_m33(){ - return m33; -} +double Metric::get_m33(){ return m33; } /*********class Node*********/ @@ -337,7 +301,8 @@ Filler::Filler(){} Filler::~Filler(){} -void Filler::treat_model(){ +void Filler::treat_model() +{ GRegion* gr; GModel* model = GModel::current(); GModel::riter it; @@ -351,8 +316,8 @@ void Filler::treat_model(){ } } -void Filler::treat_region(GRegion* gr){ - +void Filler::treat_region(GRegion* gr) +{ int NumSmooth = CTX::instance()->mesh.smoothCrossField; std::cout << "NumSmooth = " << NumSmooth << std::endl ; if(NumSmooth && (gr->dim() == 3)){ @@ -364,7 +329,6 @@ void Filler::treat_region(GRegion* gr){ Frame_field::saveCrossField("cross1.pos",scale); } -#if defined(HAVE_RTREE) unsigned int i; int j; int count; @@ -390,7 +354,7 @@ void Filler::treat_region(GRegion* gr){ std::list<GFace*>::iterator it2; std::map<MVertex*,int>::iterator it3; RTree<Node*,double,3,double> rtree; - + Frame_field::init_region(gr); Size_field::init_region(gr); Size_field::solve(gr); @@ -402,7 +366,7 @@ void Filler::treat_region(GRegion* gr){ faces.clear(); limits.clear(); - faces = gr->faces(); + faces = gr->faces(); for(it2=faces.begin();it2!=faces.end();it2++){ gf = *it2; limit = code(gf->tag()); @@ -415,7 +379,7 @@ void Filler::treat_region(GRegion* gr){ } } } - + /*for(i=0;i<gr->getNumMeshElements();i++){ element = gr->getMeshElement(i); for(j=0;j<element->getNumVertices();j++){ @@ -429,63 +393,63 @@ void Filler::treat_region(GRegion* gr){ boundary_vertices.push_back(*it); } } - + for(it=temp.begin();it!=temp.end();it++){ if((*it)->onWhat()->dim()==1){ boundary_vertices.push_back(*it); } } - + for(it=temp.begin();it!=temp.end();it++){ if((*it)->onWhat()->dim()==2){ boundary_vertices.push_back(*it); } } - + /*for(it=temp.begin();it!=temp.end();it++){ if((*it)->onWhat()->dim()<3){ boundary_vertices.push_back(*it); } }*/ //std::ofstream file("nodes.pos"); - //file << "View \"test\" {\n"; + //file << "View \"test\" {\n"; for(i=0;i<boundary_vertices.size();i++){ x = boundary_vertices[i]->x(); y = boundary_vertices[i]->y(); z = boundary_vertices[i]->z(); - + node = new Node(SPoint3(x,y,z)); compute_parameters(node,gr); node->set_layer(0); - + it3 = limits.find(boundary_vertices[i]); node->set_limit(it3->second); - + rtree.Insert(node->min,node->max,node); fifo.push(node); //print_node(node,file); } - + count = 1; while(!fifo.empty()){ parent = fifo.front(); fifo.pop(); garbage.push_back(parent); - + if(parent->get_limit()!=-1 && parent->get_layer()>=parent->get_limit()){ continue; } - + spawns.clear(); spawns.resize(6); - + for(i=0;i<6;i++){ spawns[i] = new Node(); } - + create_spawns(gr,octree,parent,spawns); - + for(i=0;i<6;i++){ ok2 = 0; individual = spawns[i]; @@ -493,18 +457,18 @@ void Filler::treat_region(GRegion* gr){ x = point.x(); y = point.y(); z = point.z(); - + if(inside_domain(octree,x,y,z)){ compute_parameters(individual,gr); individual->set_layer(parent->get_layer()+1); individual->set_limit(parent->get_limit()); - + if(far_from_boundary(octree,individual)){ wrapper.set_ok(1); wrapper.set_individual(individual); wrapper.set_parent(parent); rtree.Search(individual->min,individual->max,rtree_callback,&wrapper); - + if(wrapper.get_ok()){ fifo.push(individual); rtree.Insert(individual->min,individual->max,individual); @@ -515,16 +479,16 @@ void Filler::treat_region(GRegion* gr){ } } } - + if(!ok2) delete individual; } - + if(count%100==0){ printf("%d\n",count); } count++; } - + //file << "};\n"; int option = CTX::instance()->mesh.algo3d; @@ -538,7 +502,7 @@ void Filler::treat_region(GRegion* gr){ MeshDelaunayVolume(regions); CTX::instance()->mesh.algo3d = option; - + for(i=0;i<garbage.size();i++) delete garbage[i]; for(i=0;i<new_vertices.size();i++) delete new_vertices[i]; new_vertices.clear(); @@ -546,10 +510,10 @@ void Filler::treat_region(GRegion* gr){ rtree.RemoveAll(); Size_field::clear(); Frame_field::clear(); -#endif } -Metric Filler::get_metric(double x,double y,double z){ +Metric Filler::get_metric(double x,double y,double z) +{ Metric m; STensor3 m2; if(CTX::instance()->mesh.smoothCrossField){ @@ -573,17 +537,18 @@ Metric Filler::get_metric(double x,double y,double z){ return m; } -Metric Filler::get_metric(double x,double y,double z,GEntity* ge){ +Metric Filler::get_metric(double x,double y,double z,GEntity* ge) +{ Metric m; SMetric3 temp; SVector3 v1,v2,v3; Field* field; FieldManager* manager; - + v1 = SVector3(1.0,0.0,0.0); v2 = SVector3(0.0,1.0,0.0); v3 = SVector3(0.0,0.0,1.0); - + manager = ge->model()->getFields(); if(manager->getBackgroundField()>0){ field = manager->get(manager->getBackgroundField()); @@ -591,27 +556,29 @@ Metric Filler::get_metric(double x,double y,double z,GEntity* ge){ (*field)(x,y,z,temp,ge); } } - + m.set_m11(v1.x()); m.set_m21(v1.y()); m.set_m31(v1.z()); - + m.set_m12(v2.x()); m.set_m22(v2.y()); m.set_m32(v2.z()); - + m.set_m13(v3.x()); m.set_m23(v3.y()); m.set_m33(v3.z()); - + return m; } double Filler::get_size(double x,double y,double z){ + return Size_field::search(x,y,z); } -double Filler::get_size(double x,double y,double z,GEntity* ge){ +double Filler::get_size(double x,double y,double z,GEntity* ge) +{ double h; Field* field; FieldManager* manager; @@ -628,14 +595,16 @@ double Filler::get_size(double x,double y,double z,GEntity* ge){ return h; } -bool Filler::inside_domain(MElementOctree* octree,double x,double y,double z){ +bool Filler::inside_domain(MElementOctree* octree,double x,double y,double z) +{ MElement* element; element = (MElement*)octree->find(x,y,z,3,true); if(element!=NULL) return 1; else return 0; } -bool Filler::far_from_boundary(MElementOctree* octree,Node* node){ +bool Filler::far_from_boundary(MElementOctree* octree,Node* node) +{ double x,y,z; double h; SPoint3 point; @@ -658,7 +627,8 @@ bool Filler::far_from_boundary(MElementOctree* octree,Node* node){ else return 0; } -void Filler::compute_parameters(Node* node,GEntity* ge){ +void Filler::compute_parameters(Node* node,GEntity* ge) +{ double x,y,z; double h; Metric m; @@ -681,7 +651,9 @@ void Filler::compute_parameters(Node* node,GEntity* ge){ node->max[2] = z + sqrt3*h; } -void Filler::create_spawns(GEntity* ge,MElementOctree* octree,Node* node,std::vector<Node*>& spawns){ +void Filler::create_spawns(GEntity* ge,MElementOctree* octree, + Node* node,std::vector<Node*>& spawns) +{ double x,y,z; double x1,y1,z1; double x2,y2,z2; @@ -739,7 +711,9 @@ void Filler::create_spawns(GEntity* ge,MElementOctree* octree,Node* node,std::ve *spawns[5] = Node(SPoint3(x6,y6,z6)); } -double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,double h1,SVector3 direction){ +double Filler::improvement(GEntity* ge,MElementOctree* octree, + SPoint3 point,double h1,SVector3 direction) +{ double x,y,z; double average; double h2; @@ -748,7 +722,7 @@ double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,doub x = point.x() + h1*direction.x(); y = point.y() + h1*direction.y(); z = point.z() + h1*direction.z(); - + if(inside_domain(octree,x,y,z)){ h2 = get_size(x,y,z); } @@ -756,14 +730,14 @@ double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,doub coeffA = 1.0; coeffB = 0.16; - + if(h2>h1){ average = coeffA*h1 + (1.0-coeffA)*h2; } else{ average = coeffB*h1 + (1.0-coeffB)*h2; } - + return average; } @@ -771,11 +745,11 @@ int Filler::code(int tag){ int limit; std::string s; std::stringstream s2; - + limit = -1; s2 << tag; s = s2.str(); - + if(s.length()>=5){ if(s.at(0)=='1' && s.at(1)=='1' && s.at(2)=='1' && s.at(3)=='1' && s.at(4)=='1'){ limit = 0; @@ -784,26 +758,30 @@ int Filler::code(int tag){ limit = 1; } } - + return limit; } -int Filler::get_nbr_new_vertices(){ +int Filler::get_nbr_new_vertices() +{ return new_vertices.size(); } -MVertex* Filler::get_new_vertex(int i){ +MVertex* Filler::get_new_vertex(int i) +{ return new_vertices[i]; } -void Filler::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){ +void Filler::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file) +{ file << "SL (" << p1.x() << ", " << p1.y() << ", " << p1.z() << ", " << p2.x() << ", " << p2.y() << ", " << p2.z() << ")" << "{10, 20};\n"; } -void Filler::print_node(Node* node,std::ofstream& file){ +void Filler::print_node(Node* node,std::ofstream& file) +{ double x,y,z; double x1,y1,z1; double x2,y2,z2; diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp index 9d4361d4143fb36ecfc5cdb445a11c54457a083a..bb88b2519101cd1863735edd27d629fda08720b0 100644 --- a/Mesh/surfaceFiller.cpp +++ b/Mesh/surfaceFiller.cpp @@ -13,20 +13,15 @@ #include "OS.h" #include <queue> #include <stack> - -/// Here, we aim at producing a set of points that -/// enables to generate a nice quad mesh - -#if defined(HAVE_RTREE) #include "rtree.h" -#endif - #include "MVertex.h" #include "MElement.h" -//#include "directions3D.h" #include "BackgroundMesh.h" #include "intersectCurveSurface.h" +// Here, we aim at producing a set of points that enables to generate a nice +// quad mesh + using namespace std; static const double FACTOR = .71; @@ -37,15 +32,10 @@ static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.}; //static const double DIRS [NUMDIR] = {0.0}; // END PE MODIF - - -/// a rectangle in the tangent plane is transformed -/// into a parallelogram. We define an exclusion zone -/// that is centered around a vertex and that is used -/// in a r-tree structure for generating points with the -/// right spacing in the tangent plane - -#if defined(HAVE_RTREE) +// a rectangle in the tangent plane is transformed into a parallelogram. We +// define an exclusion zone that is centered around a vertex and that is used in +// a r-tree structure for generating points with the right spacing in the +// tangent plane struct surfacePointWithExclusionRegion { MVertex *_v; @@ -62,9 +52,11 @@ struct surfacePointWithExclusionRegion { | + p1 -*/ - - surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], SPoint2 &_mp, SMetric3 & meshMetric, surfacePointWithExclusionRegion *father = 0){ + */ + surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], + SPoint2 &_mp, SMetric3 & meshMetric, + surfacePointWithExclusionRegion *father = 0) + { _v = v; _meshMetric = meshMetric; _center = _mp; @@ -83,10 +75,11 @@ struct surfacePointWithExclusionRegion { } //test !!! // const double s = backgroundMesh::current()->getSmoothness(_mp.x(),_mp.y(),0); - // if(s > 0.02) + // if(s > 0.02) // _distanceSummed = 10000*s; } - bool inExclusionZone (const SPoint2 &p){ + bool inExclusionZone (const SPoint2 &p) + { double mat[2][2]; double b[2] , uv[2]; mat[0][0]= _q[1].x()-_q[0].x(); @@ -109,13 +102,15 @@ struct surfacePointWithExclusionRegion { if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true; return false; } - void minmax (double _min[2], double _max[2]) const{ + void minmax (double _min[2], double _max[2]) const + { _min[0] = std::min(std::min(std::min(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x()); _min[1] = std::min(std::min(std::min(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y()); _max[0] = std::max(std::max(std::max(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x()); _max[1] = std::max(std::max(std::max(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y()); } - void print (FILE *f, int i){ + void print (FILE *f, int i) + { fprintf(f,"SP(%g,%g,%g){%d};\n",_center.x(),_center.y(),0.0,i); fprintf(f,"SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n", _q[0].x(),_q[0].y(),0.0, @@ -140,7 +135,8 @@ struct smoothness_point_pair{ class compareSurfacePointWithExclusionRegionPtr_Smoothness { public: - inline bool operator () (const smoothness_point_pair &a, const smoothness_point_pair &b) const + inline bool operator () (const smoothness_point_pair &a, + const smoothness_point_pair &b) const { if (a.rank == b.rank){ if(a.ptr->_distanceSummed > b.ptr->_distanceSummed) return false; @@ -152,11 +148,11 @@ class compareSurfacePointWithExclusionRegionPtr_Smoothness } }; - class compareSurfacePointWithExclusionRegionPtr { public: - inline bool operator () (const surfacePointWithExclusionRegion *a, const surfacePointWithExclusionRegion *b) const + inline bool operator () (const surfacePointWithExclusionRegion *a, + const surfacePointWithExclusionRegion *b) const { if(a->_distanceSummed > b->_distanceSummed) return false; if(a->_distanceSummed < b->_distanceSummed) return true; @@ -164,10 +160,8 @@ class compareSurfacePointWithExclusionRegionPtr } }; - - - -bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){ +bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point) +{ my_wrapper *w = static_cast<my_wrapper*>(point); if (neighbour->inExclusionZone(w->_p)){ @@ -180,7 +174,8 @@ bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){ bool inExclusionZone (SPoint2 &p, RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree, - std::vector<surfacePointWithExclusionRegion*> & all ){ + std::vector<surfacePointWithExclusionRegion*> & all) +{ // should assert that the point is inside the domain if (!backgroundMesh::current()->inDomain(p.x(),p.y(),0)) return true; @@ -192,16 +187,14 @@ bool inExclusionZone (SPoint2 &p, for (unsigned int i=0;i<all.size();++i){ if (all[i]->inExclusionZone(p)){ - // printf("%g %g is in exclusion zone of %g %g\n",p.x(),p.y(),all[i]._center.x(),all[i]._center.y()); + // printf("%g %g is in exclusion zone of %g %g\n", + // p.x(),p.y(),all[i]._center.x(),all[i]._center.y()); return true; } } return false; } - - - // assume a point on the surface, compute the 4 possible neighbors. // // ^ t2 @@ -224,7 +217,6 @@ bool compute4neighbors (GFace *gf, // the surface SPoint2 newP[4][NUMDIR], // look into other directions SMetric3 &metricField, FILE *crossf = 0) // the mesh metric { - //Range<double> rangeU = gf->parBounds(0); //Range<double> rangeV = gf->parBounds(1); @@ -248,7 +240,6 @@ bool compute4neighbors (GFace *gf, // the surface } } - // get the unit normal at that point Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(midpoint[0],midpoint[1])); SVector3 s1 = der.first(); @@ -260,7 +251,6 @@ bool compute4neighbors (GFace *gf, // the surface double N = dot(s2,s2); double E = dot(s1,s2); - // compute the first fundamental form i.e. the metric tensor at the point // M_{ij} = s_i \cdot s_j double metric[2][2] = {{M,E},{E,N}}; @@ -271,7 +261,8 @@ bool compute4neighbors (GFace *gf, // the surface SVector3 basis_v = crossprod(n,basis_u); for (int DIR = 0 ; DIR < NUMDIR ; DIR ++){ - double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) + DIRS[DIR]; + double quadAngle = backgroundMesh::current()->getAngle + (midpoint[0], midpoint[1],0) + DIRS[DIR]; // normalize vector t1 that is tangent to gf at midpoint SVector3 t1 = basis_u * cos (quadAngle) + basis_v * sin (quadAngle) ; @@ -280,10 +271,18 @@ bool compute4neighbors (GFace *gf, // the surface // compute the second direction t2 and normalize (t1,t2,n) is the tangent frame SVector3 t2 = crossprod(n,t1); t2.normalize(); - if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),t1.x(),t1.y(),t1.z()); - if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),t2.x(),t2.y(),t2.z()); - if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),-t1.x(),-t1.y(),-t1.z()); - if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),-t2.x(),-t2.y(),-t2.z()); + if (DIR == 0 && crossf) + fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n", + v_center->x(),v_center->y(),v_center->z(),t1.x(),t1.y(),t1.z()); + if (DIR == 0 && crossf) + fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n", + v_center->x(),v_center->y(),v_center->z(),t2.x(),t2.y(),t2.z()); + if (DIR == 0 && crossf) + fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n", + v_center->x(),v_center->y(),v_center->z(),-t1.x(),-t1.y(),-t1.z()); + if (DIR == 0 && crossf) + fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n", + v_center->x(),v_center->y(),v_center->z(),-t2.x(),-t2.y(),-t2.z()); double size_1 = sqrt(1. / dot(t1,metricField,t1)); double size_2 = sqrt(1. / dot(t2,metricField,t2)); @@ -296,21 +295,22 @@ bool compute4neighbors (GFace *gf, // the surface double rhs1[2] = {dot(t1,s1),dot(t1,s2)}, covar1[2]; bool singular = false; if (!sys2x2(metric,rhs1,covar1)){ - Msg::Info("Argh surface %d %g %g %g -- %g %g %g -- %g %g",gf->tag(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z(),size_1,size_2); + Msg::Info("Argh surface %d %g %g %g -- %g %g %g -- %g %g", + gf->tag(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z(),size_1,size_2); covar1[1] = 1.0; covar1[0] = 0.0; singular = true; } double rhs2[2] = {dot(t2,s1),dot(t2,s2)}, covar2[2]; if (!sys2x2(metric,rhs2,covar2)){ - Msg::Info("Argh surface %d %g %g %g -- %g %g %g",gf->tag(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z()); + Msg::Info("Argh surface %d %g %g %g -- %g %g %g", + gf->tag(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z()); covar2[0] = 1.0; covar2[1] = 0.0; singular = true; } - // transform the sizes with respect to the metric - // consider a vector v of size 1 in the parameter plane - // its length is sqrt (v^T M v) --> if I want a real size - // of size1 in direction v, it should be sqrt(v^T M v) * size1 + // transform the sizes with respect to the metric consider a vector v of + // size 1 in the parameter plane its length is sqrt (v^T M v) --> if I want + // a real size of size1 in direction v, it should be sqrt(v^T M v) * size1 double l1 = sqrt(covar1[0]*covar1[0]+covar1[1]*covar1[1]); double l2 = sqrt(covar2[0]*covar2[0]+covar2[1]*covar2[1]); @@ -378,17 +378,17 @@ bool compute4neighbors (GFace *gf, // the surface // printf("L = %12.5E D = %12.5E ERR = %12.5E\n",L,D,100*fabs(D-L)/(D+L)); } - if (1 && goNonLinear){//---------------------------------------------------// - surfaceFunctorGFace ss (gf); // - SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1*(1.0),t2*(1.0)}; // - for (int i=0;i<4;i++){ // + if (1 && goNonLinear){ + surfaceFunctorGFace ss (gf); + SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1*(1.0),t2*(1.0)}; + for (int i=0;i<4;i++){ if (ERR[i] > 12){ - double uvt[3] = {newPoint[i][0],newPoint[i][1],0.0}; // + double uvt[3] = {newPoint[i][0],newPoint[i][1],0.0}; // printf("Intersecting with circle N = %g %g %g dir = %g %g %g R = %g p = %g %g %g\n",n.x(),n.y(),n.z(),dirs[i].x(),dirs[i].y(),dirs[i].z(),L,v_center->x(),v_center->y(),v_center->z()); curveFunctorCircle cf (dirs[i],n, SVector3(v_center->x(),v_center->y(),v_center->z()), L); - if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-3)){ // + if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-3)){ GPoint pp = gf->point(SPoint2(uvt[0],uvt[1])); double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) + (pp.y() - v_center->y())*(pp.y() - v_center->y()) + @@ -396,15 +396,15 @@ bool compute4neighbors (GFace *gf, // the surface double DP = sqrt ((newPoint[i][0]-uvt[0])*(newPoint[i][0]-uvt[0]) + (newPoint[i][1]-uvt[1])*(newPoint[i][1]-uvt[1])); double newErr = 100*fabs(D-L)/(D+L); - // if (v_center->onWhat() != gf && gf->tag() == 3){ - // crossField2d::normalizeAngle (uvt[2]); - // printf("INTERSECT angle = %g DP %g\n",uvt[2],DP); - // } + // if (v_center->onWhat() != gf && gf->tag() == 3){ + // crossField2d::normalizeAngle (uvt[2]); + // printf("INTERSECT angle = %g DP %g\n",uvt[2],DP); + // } if (newErr < 1 && DP < .1){ // printf("%12.5E vs %12.5E : %12.5E %12.5E vs %12.5E %12.5E \n",ERR[i],newErr,newPoint[i][0],newPoint[i][1],uvt[0],uvt[1]); - newPoint[i][0] = uvt[0]; // - newPoint[i][1] = uvt[1]; // - } // + newPoint[i][0] = uvt[0]; + newPoint[i][1] = uvt[1]; + } // printf("OK\n"); } else{ @@ -412,8 +412,8 @@ bool compute4neighbors (GFace *gf, // the surface // printf("NOT OK\n"); } } - } // - } /// end non linear -------------------------------------------------// + } + } // end non linear // return the four new vertices for (int i=0;i<4;i++){ @@ -423,10 +423,10 @@ bool compute4neighbors (GFace *gf, // the surface return true; } -// --------------------------------------------------------------------------------------------- - // recover element around vertex v and interpolate smoothness on this element... -double get_smoothness(MVertex *v, GFace *gf, const map<MVertex*,double> &vertices2smoothness){ +double get_smoothness(MVertex *v, GFace *gf, + const map<MVertex*,double> &vertices2smoothness) +{ // recover element around MVertex v //cout << "Looking for element around point (" << v->x() << "," << v->y() << "," << v->z() << ")" << endl; SPoint3 sp3(v->x(), v->y(), v->z()); @@ -475,9 +475,8 @@ double get_smoothness(MVertex *v, GFace *gf, const map<MVertex*,double> &vertice return res; } -// --------------------------------------------------------------------------------------------- - -void print_nodal_info_int(string filename, map<MVertex*, int> &mapp){ +void print_nodal_info_int(string filename, map<MVertex*, int> &mapp) +{ ofstream out(filename.c_str()); out << "View \"\"{" << endl; @@ -490,9 +489,8 @@ void print_nodal_info_int(string filename, map<MVertex*, int> &mapp){ out.close(); } -// --------------------------------------------------------------------------------------------- - -void print_nodal_info_double(string filename, map<MVertex*, double> &mapp){ +void print_nodal_info_double(string filename, map<MVertex*, double> &mapp) +{ ofstream out(filename.c_str()); out << "View \"\"{" << endl; @@ -505,9 +503,8 @@ void print_nodal_info_double(string filename, map<MVertex*, double> &mapp){ out.close(); } -// --------------------------------------------------------------------------------------------- - -void export_point(surfacePointWithExclusionRegion *sp, int DIR, FILE *crossf, GFace *gf){ +void export_point(surfacePointWithExclusionRegion *sp, int DIR, FILE *crossf, GFace *gf) +{ // get the unit normal at that point Pair<SVector3, SVector3> der = gf->firstDer(sp->_center); SVector3 s1 = der.first(); @@ -555,9 +552,12 @@ void export_point(surfacePointWithExclusionRegion *sp, int DIR, FILE *crossf, GF fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",sp->_v->x(),sp->_v->y(),sp->_v->z(),-t2.x()*size_2,-t2.y()*size_2,-t2.z()*size_2); } -// --------------------------------------------------------------------------------------------- - -bool get_local_sizes_and_directions(const MVertex *v_center, const SPoint2 &midpoint, const int DIR, GFace* gf, double (&covar1)[2], double (&covar2)[2], double &size_param_1, double &size_param_2, double &L, SVector3 &t1, SVector3 &t2, SVector3 &n, FILE *crossf=NULL){ +bool get_local_sizes_and_directions(const MVertex *v_center, const SPoint2 &midpoint, + const int DIR, GFace* gf, double (&covar1)[2], + double (&covar2)[2], double &size_param_1, + double &size_param_2, double &L, SVector3 &t1, + SVector3 &t2, SVector3 &n, FILE *crossf=NULL) +{ //bool get_RK_stuff(const MVertex *v_center, const SPoint2 &midpoint, const int DIR, GFace* gf, double (&covar1)[2], double (&covar2)[2], double &size_param_1, double &size_param_2, double &L, SVector3 &t1, SVector3 &t2, SVector3 &n, FILE *crossf, const SVector3 &previous_t1, const SVector3 &previous_t2, bool use_previous_basis=false, bool export_cross=true){ // !!!!!!!!!!!! check if point is in domain (for RK !!!) @@ -692,16 +692,12 @@ bool get_local_sizes_and_directions(const MVertex *v_center, const SPoint2 &midp return true; } -#endif - -// --------------------------------------------------------------------------------------------- - // using fifo based on smoothness criteria -void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics){ +void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed, + std::vector<SMetric3> &metrics) +{ cout << endl << "------------------------------------------" << endl << " PACKINGOFPARALLELOGRAMS: NEW ALGO BASED ON SMOOTHNESS" << endl << "------------------------------------------" << endl; -#if defined(HAVE_RTREE) const bool goNonLinear = true; - const bool debug = false; // build vertex -> neighbors table @@ -913,21 +909,17 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed } fprintf(f,"};"); fclose(f); -#endif } - -// --------------------------------------------------------------------------------------------- - - // fills a surface with points in order to build a nice // quad mesh ------------ -void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics){ +void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, + std::vector<SMetric3> &metrics) +{ //PE MODIF // packingOfParallelogramsSmoothness(gf,packed,metrics); // return; // END PE MODIF -#if defined(HAVE_RTREE) const bool goNonLinear = true; @@ -1038,5 +1030,4 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec fclose(f); // printf("packed.size = %d\n",packed.size()); // delete rtree; -#endif } diff --git a/contrib/rtree/README.TXT b/contrib/rtree/README.TXT deleted file mode 100644 index a2f23a279f25855d60587e27dc9d51b3a7726c8a..0000000000000000000000000000000000000000 --- a/contrib/rtree/README.TXT +++ /dev/null @@ -1,43 +0,0 @@ - -TITLE - - R-TREES: A DYNAMIC INDEX STRUCTURE FOR SPATIAL SEARCHING - -DESCRIPTION - - A C++ templated version of the RTree algorithm. - For more information please read the comments in RTree.h - -AUTHORS - - * 1983 Original algorithm and test code by Antonin Guttman and Michael Stonebraker, UC Berkely - * 1994 ANCI C ported from original test code by Melinda Green - melinda@superliminal.com - * 1995 Sphere volume fix for degeneracy problem submitted by Paul Brook - * 2004 Templated C++ port by Greg Douglas - -LICENSE: - - Entirely free for all uses. Enjoy! - -FILES - * RTree.h The C++ templated RTree implementation. Well commented. - * Test.cpp A simple test program, ported from the original C version. - * MemoryTest.cpp A more rigourous test to validate memory use. - * README.TXT This file. - -TO BUILD - - To build a test, compile only one of the test files with RTree.h. - Both test files contain a main() function. - -RECENT CHANGE LOG - -05 Jan 2010 -o Fixed Iterator GetFirst() - Previous fix was not incomplete - -03 Dec 2009 -o Added Iteartor GetBounds() -o Added Iterator usage to simple test -o Fixed Iterator GetFirst() - Thanks Mathew Riek -o Minor updates for MSVC 2005/08 compilers - diff --git a/doc/texinfo/cmake_options.texi b/doc/texinfo/cmake_options.texi index 745c6207848af3218442e5acab7e7a7d5a21abaa..c057a6667f5f3515cf1ca83ac968997aef15f22c 100644 --- a/doc/texinfo/cmake_options.texi +++ b/doc/texinfo/cmake_options.texi @@ -29,6 +29,8 @@ Enable CGNS mesh export (experimental) (default: OFF) Enable Cairo to render fonts (experimental) (default: ON) @item ENABLE_CHACO Enable Chaco mesh partitioner (alternative to Metis) (default: ON) +@item ENABLE_COMPRESSED_IO +Enable compressed (gzip) input/output using zlib (default: OFF) @item ENABLE_DINTEGRATION Enable discrete integration (needed for levelsets) (default: ON) @item ENABLE_FLTK @@ -57,16 +59,24 @@ Enable built-in MPEG movie encoder (default: ON) Enable MPI (mostly for parser and solver - mesh generation is sequential) (default: OFF) @item ENABLE_MSVC_STATIC_RUNTIME Enable static Visual C++ runtime (default: OFF) +@item ENABLE_MUMPS +Enable MUMPS sparse direct linear solver (default: OFF) @item ENABLE_NATIVE_FILE_CHOOSER Enable native file chooser in GUI (default: ON) @item ENABLE_NETGEN Enable Netgen 3D frontal mesh generator (default: ON) +@item ENABLE_NUMPY +Enable conversion between fullMatrix and numpy array (requires SWIG) (default: OFF) +@item ENABLE_PETSC4PY +Enable petsc4py wrappers for petsc matrices (default: ON) @item ENABLE_OCC Enable Open CASCADE geometrical models (default: ON) @item ENABLE_ONELAB -Enable OneLab solver interface (default: ON) +Enable ONELAB solver interface (default: ON) +@item ENABLE_ONELAB2 +Enable experimental ONELAB-Cloud solver interface (default: OFF) @item ENABLE_ONELAB_METAMODEL -Enable OneLab metamodels (experimental) (default: ON) +Enable ONELAB metamodels (experimental) (default: ON) @item ENABLE_OPENMP Enable OpenMP (experimental) (default: OFF) @item ENABLE_OPTHOM @@ -99,6 +109,8 @@ Enable SLEPc eigensolvers (required for conformal compounds) (default: ON) Enable built-in finite element solvers (required for compounds) (default: ON) @item ENABLE_TAUCS Enable Taucs linear solver (alternative to PETSc) (default: ON) +@item ENABLE_TCMALLOC +Enable libtcmalloc, a fast malloc implementation but that does not release memory (default: OFF) @item ENABLE_TETGEN Enable Tetgen 3D initial mesh generator (default: ON) @item ENABLE_TETGEN_OLD diff --git a/doc/texinfo/opt_plugin.texi b/doc/texinfo/opt_plugin.texi index 1c2f999605c8064f56c82315e80453a413f3b36b..08c2b0dd090d4f8720cce97bb1eb17ef4648e665 100644 --- a/doc/texinfo/opt_plugin.texi +++ b/doc/texinfo/opt_plugin.texi @@ -724,6 +724,20 @@ Default value: @code{0} Default value: @code{-1} @end table +@item Plugin(MeshSubEntities) +Plugin(MeshSubEntities) creates mesh elements for the entities of dimension `OutputDimension' (0 for vertices, 1 for edges, 2 for faces) of the `InputPhysicalGroup' of dimension `InputDimension'. The plugin creates new elements belonging to `OutputPhysicalGroup'. +Numeric options: +@table @code +@item InputDimension +Default value: @code{1} +@item InputPhysicalGroup +Default value: @code{1} +@item OuputDimension +Default value: @code{0} +@item OuputPhysicalGroup +Default value: @code{2000} +@end table + @item Plugin(MinMax) Plugin(MinMax) computes the min/max of a view.@* @*