From bf2cdc51c7aaad44929446d7d1b2e4ef7c8c184b Mon Sep 17 00:00:00 2001 From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be> Date: Wed, 31 Oct 2001 16:33:46 +0000 Subject: [PATCH] *** empty log message *** --- DataStr/Tree.cpp | 9 +- DataStr/Tree.h | 9 +- DataStr/avl.cpp | 31 +--- DataStr/avl.h | 29 ++- Fltk/Socket.cpp | 5 +- Mesh/3D_Mesh.cpp | 91 ++++++++-- Mesh/mGeomSearch.h | 184 ++++++++++++++++++++ Plugin/Makefile | 4 +- Plugin/Plugin.cpp | 3 +- Triangle/Makefile | 4 +- benchmarks/3d/magnetron/magnetron1-full.geo | 2 +- 11 files changed, 313 insertions(+), 58 deletions(-) create mode 100644 Mesh/mGeomSearch.h diff --git a/DataStr/Tree.cpp b/DataStr/Tree.cpp index 091c602077..a154f3f06f 100644 --- a/DataStr/Tree.cpp +++ b/DataStr/Tree.cpp @@ -1,4 +1,4 @@ -// $Id: Tree.cpp,v 1.7 2001-04-08 20:36:49 geuzaine Exp $ +// $Id: Tree.cpp,v 1.8 2001-10-31 16:33:46 remacle Exp $ #include <stdlib.h> #include <string.h> @@ -167,13 +167,6 @@ int Tree_Right(Tree_T *tree, void *data) return (1); } -void Tree_Action(Tree_T *tree, void (*action) (void *data, void *dummy)) -{ - if(!tree) return; - - avl_foreach(tree->root, action, AVL_FORWARD); -} - int Tree_Size(Tree_T *tree) { if(!tree) return 0; diff --git a/DataStr/Tree.h b/DataStr/Tree.h index 85be99ed56..31f6f158ff 100644 --- a/DataStr/Tree.h +++ b/DataStr/Tree.h @@ -21,8 +21,15 @@ void *Tree_PQuery(Tree_T *Tree, void *data); int Tree_Suppress(Tree_T *Tree, void *data); int Tree_Left(Tree_T *tree, void *data); int Tree_Right(Tree_T *tree, void *data); -void Tree_Action(Tree_T *tree, void (*action) (void *data, void *dummy)); int Tree_Size(Tree_T *tree) ; +inline void Tree_Action(Tree_T *tree, void (*action) (void *data, void *dummy)) +{ + if(!tree) return; + + avl_foreach(tree->root, action, AVL_FORWARD); +} + + #endif diff --git a/DataStr/avl.cpp b/DataStr/avl.cpp index c10d4b2ce9..f30b67579e 100644 --- a/DataStr/avl.cpp +++ b/DataStr/avl.cpp @@ -1,4 +1,4 @@ -// $Id: avl.cpp,v 1.5 2001-01-08 08:05:41 geuzaine Exp $ +// $Id: avl.cpp,v 1.6 2001-10-31 16:33:46 remacle Exp $ /* * This is a modified version for Gmsh (for c++, 64-bit architectures, etc.) @@ -35,7 +35,6 @@ #include "Malloc.h" #define ALLOC(type, number) (type *) Malloc((unsigned) sizeof(type) * number) -#define NIL(type) (type *) 0 #define FREE(item) (void) Free(item) #define XRNMAX(a,b) ((a) > (b) ? (a) : (b)) #define HEIGHT(node) (node == NIL(avl_node) ? -1 : (node)->height) @@ -57,8 +56,6 @@ static avl_node *find_rightmost(avl_node **node_p); static void do_rebalance(avl_node ***stack_nodep, int stack_n); static void rotate_left(avl_node **node_p); static void rotate_right(avl_node **node_p); -static void avl_walk_forward(avl_node *node, void (*func)(void *key, void *value)); -static void avl_walk_backward(avl_node *node, void (*func)(void *key, void *value)); static void free_entry(avl_node *node, void (*key_free)(void *key), void (*value_free)(void *value)); static avl_node *new_node(void *key, void *value); @@ -307,32 +304,6 @@ static void rotate_right(avl_node **node_p) compute_height(new_root); } -static void avl_walk_forward(avl_node *node, void (*func)(void *key, void *value)) -{ - if (node != NIL(avl_node)) { - avl_walk_forward(node->left, func); - (*func)(node->key, node->value); - avl_walk_forward(node->right, func); - } -} - -static void avl_walk_backward(avl_node *node, void (*func)(void *key, void *value)) -{ - if (node != NIL(avl_node)) { - avl_walk_backward(node->right, func); - (*func)(node->key, node->value); - avl_walk_backward(node->left, func); - } -} - -void avl_foreach(avl_tree *tree, void (*func)(void *key, void *value), int direction) -{ - if (direction == AVL_FORWARD) { - avl_walk_forward(tree->root, func); - } else { - avl_walk_backward(tree->root, func); - } -} int avl_extremum(avl_tree *tree, int side, void **value_p) { diff --git a/DataStr/avl.h b/DataStr/avl.h index 6ff601f6bb..f343449f12 100644 --- a/DataStr/avl.h +++ b/DataStr/avl.h @@ -59,17 +59,44 @@ struct avl_generator_struct { #define AVL_MOST_RIGHT 1 #define avl_is_member(tree, key) avl_lookup(tree, key, (void **) 0) +#define NIL(type) (type *) 0 #define avl_foreach_item(table, gen, dir, key_p, value_p) \ for(gen = avl_init_gen(table, dir); \ avl_gen(gen, key_p, value_p) || (avl_free_gen(gen),0);) +inline void avl_walk_forward(avl_node *node, void (*func)(void *key, void *value)) +{ + if (node != NIL(avl_node)) { + avl_walk_forward(node->left, func); + (*func)(node->key, node->value); + avl_walk_forward(node->right, func); + } +} + +inline void avl_walk_backward(avl_node *node, void (*func)(void *key, void *value)) +{ + if (node != NIL(avl_node)) { + avl_walk_backward(node->right, func); + (*func)(node->key, node->value); + avl_walk_backward(node->left, func); + } +} + +inline void avl_foreach(avl_tree *tree, void (*func)(void *key, void *value), int direction) +{ + if (direction == AVL_FORWARD) { + avl_walk_forward(tree->root, func); + } else { + avl_walk_backward(tree->root, func); + } +} + avl_tree *avl_init_table(int (*compar)(const void *key1, const void *key2)); int avl_lookup(avl_tree *tree, void *key, void **value_p); int avl_insert(avl_tree *tree, void *key, void *value); int avl_delete(avl_tree *tree, void **key_p, void **value_p); -void avl_foreach(avl_tree *tree, void (*func)(void *key, void *value), int direction); void avl_free_table(avl_tree *tree, void (*key_free)(void *key), void (*value_free)(void *value)); int avl_count(avl_tree *tree); int avl_check_tree(avl_tree *tree); diff --git a/Fltk/Socket.cpp b/Fltk/Socket.cpp index f1bd75856e..b30cb5232b 100644 --- a/Fltk/Socket.cpp +++ b/Fltk/Socket.cpp @@ -1,4 +1,4 @@ -/* $Id: Socket.cpp,v 1.12 2001-08-26 12:12:18 geuzaine Exp $ */ +/* $Id: Socket.cpp,v 1.13 2001-10-31 16:33:46 remacle Exp $ */ #include <stdio.h> #include <stdlib.h> @@ -147,7 +147,8 @@ int Socket_StartProgram(char *progname, char *sockname){ } len = sizeof(from); - if ( (sock = accept(s, (struct sockaddr *)&from, &len)) < 0) { + socklen_t slen = len; + if ( (sock = accept(s, (struct sockaddr *)&from, &slen)) < 0) { Msg(GERROR, "Socket accept failed"); return -1; } diff --git a/Mesh/3D_Mesh.cpp b/Mesh/3D_Mesh.cpp index 86d642bbbb..d9ee2ccf48 100644 --- a/Mesh/3D_Mesh.cpp +++ b/Mesh/3D_Mesh.cpp @@ -1,4 +1,4 @@ -// $Id: 3D_Mesh.cpp,v 1.31 2001-10-29 08:52:20 geuzaine Exp $ +// $Id: 3D_Mesh.cpp,v 1.32 2001-10-31 16:33:46 remacle Exp $ /* @@ -16,12 +16,15 @@ */ +#include <time.h> +#include <vector> #include "Gmsh.h" #include "Numeric.h" #include "Mesh.h" #include "3D_Mesh.h" #include "Create.h" #include "Context.h" +#include "mGeomSearch.h" extern Mesh *THEM, *LOCAL; extern Context_T CTX; @@ -39,6 +42,7 @@ static int ZONEELIMINEE, Methode = 0; Simplex MyNewBoundary; int Alerte_Point_Scabreux; + void DebugSimplexe (Simplex * s){ int i; @@ -146,7 +150,7 @@ int Pt_In_Volume (double X, double Y, double Z, Mesh * m, return (0); } -int Pt_In_Circum (Simplex * s, Vertex * v){ +inline int Pt_In_Circum (Simplex * s, Vertex * v){ double d1, d2, eps; /* Determine si un point est dans le cercle circonscrit a un simplexe */ @@ -168,6 +172,35 @@ int Pt_In_Circum (Simplex * s, Vertex * v){ return (0); } +struct SimplexInteriorCheck +{ + bool operator() ( Simplex * s , double x[3], double u[3]) + { + Vertex v(x[0],x[1],x[2]); + return Pt_In_Circum (s, &v); + } +}; + +struct SimplexInBox +{ + double sizeBox; + SimplexInBox(double sb) : sizeBox(sb) {} + bool operator() ( Simplex * s , double min[3], double max[3]) + { + double ss; + if(sizeBox > s->Radius) ss = s->Radius; + else ss = sizeBox; + min[0] = s->Center.X - ss; + max[0] = s->Center.X + ss; + min[1] = s->Center.Y - ss; + max[1] = s->Center.Y + ss; + min[2] = s->Center.Z - ss; + max[2] = s->Center.Z + ss; + } +}; + +AOMD::mGeomSearch < Simplex, SimplexInBox, SimplexInteriorCheck > *fastSearch = 0; + void Action_First_Simplexes (void *a, void *b){ Simplex **q; @@ -467,6 +500,8 @@ int recur_bowyer (Simplex * s){ return 1; } + + bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){ int i; Simplex *s; @@ -536,7 +571,10 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){ else{ Tree_Add (m->Vertices, &THEV); for (i = 0; i < List_Nbr (Simplexes_New); i++){ - Tree_Add (m->Simplexes, List_Pointer (Simplexes_New, i)); + Simplex *theNewS; + List_Read (Simplexes_New, i, &theNewS); + Tree_Add (m->Simplexes, &theNewS); + if(fastSearch)fastSearch->AddObject(theNewS); } /* Suppression des simplexes elimines */ @@ -545,6 +583,7 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){ List_Read (Simplexes_Destroyed, i, &s); if (!Tree_Suppress (m->Simplexes, &s)) Msg(GERROR, "Impossible to delete simplex"); + if(fastSearch)fastSearch->RemoveObject(s); Free_Simplex (&s,0); } @@ -562,22 +601,50 @@ bool Bowyer_Watson (Mesh * m, Vertex * v, Simplex * S, int force){ void Convex_Hull_Mesh (List_T * Points, Mesh * m){ int i, j, N, n; int Nbr_OK = 0, Nbr_NOTOK = 0; + double xxx[3],uuu[3]; N = List_Nbr (Points); n = IMAX (N / 20, 1); + clock_t t1 = clock(); + Box_6_Tetraedron (Points, m); - // List_Sort (Points, comparePosition); + // List_Sort (Points, comparePosition); + + SimplexInteriorCheck ick; + SimplexInBox ibx (LC3D/30.); + + fastSearch = new AOMD::mGeomSearch < Simplex, SimplexInBox, SimplexInteriorCheck > + ( m->Grid.min.X,m->Grid.max.X, + m->Grid.min.Y,m->Grid.max.Y, + m->Grid.min.Z,m->Grid.max.Z, + ibx,ick, 30,30,30); + + int nbEasy = 0, nbFast = 0; for (i = 0; i < N; i++){ THES = NULL; List_Read (Points, i, &THEV); - - if (Simplexes_New) - for (j = 0; j < List_Nbr (Simplexes_New); j++){ - Action_First_Simplexes (List_Pointer (Simplexes_New, j), NULL); + xxx[0]= THEV->Pos.X; + xxx[1]= THEV->Pos.Y; + xxx[2]= THEV->Pos.Z; + if(fastSearch) THES = fastSearch->find(xxx,uuu); + if(!THES) + { + if (Simplexes_New) + for (j = 0; j < List_Nbr (Simplexes_New); j++){ + Action_First_Simplexes (List_Pointer (Simplexes_New, j), NULL); + if(THES) + { + nbEasy++; + break; + } + } + } + else + { + nbFast++; } - if (!THES){ Tree_Action (m->Simplexes, Action_First_Simplexes); Nbr_OK++; @@ -589,7 +656,7 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){ volume = 0.0; Tree_Action (m->Simplexes, VSIM); Msg(STATUS3, "Nod=%d/%d Elm=%d", i+1,N,Tree_Nbr(m->Simplexes)); - Msg(STATUS1, "Vol=%g",volume); + Msg(STATUS1, "Vol=%g (%d %d %d)",volume,nbFast,nbEasy,i-nbFast-nbEasy); } if (!THES){ Msg(WARNING, "Vertex (%g,%g,%g) in no simplex", THEV->Pos.X,THEV->Pos.Y,THEV->Pos.Z); @@ -622,6 +689,10 @@ void Convex_Hull_Mesh (List_T * Points, Mesh * m){ } } } + clock_t t2 = clock(); + + if(fastSearch)delete fastSearch; + printf("nbEasy = %d nbFast = %d N = %d t = %lf\n",nbEasy,nbFast,N,(double)(t2-t1)/CLOCKS_PER_SEC); } void suppress_vertex (void *data, void *dum){ diff --git a/Mesh/mGeomSearch.h b/Mesh/mGeomSearch.h new file mode 100644 index 0000000000..3354e8eb4b --- /dev/null +++ b/Mesh/mGeomSearch.h @@ -0,0 +1,184 @@ +#ifndef _MGEOM_SEARCH_H +#define _MGEOM_SEARCH_H +#include <vector> +#include <algorithm> +#define MIN(x,y) ((x<y)?(x):(y)) +#define MAX(x,y) ((x>y)?(x):(y)) +#define TOL 1.e-06 + +namespace AOMD { + + template <class T> + class Brick + { + public: + std::vector<T*> Objects; + Brick(){} + T* operator [] (int i) + { + if(i < 0 || i >= Objects.size())throw i; + T *obj = Objects[i]; + return obj; + } + int size(){return Objects.size();} + }; + + template <class T, class GetBox, class InteriorCheck> + class mGeomSearch + { + Brick<T> *theTable; + int Nx,Ny,Nz; + double Xmin,Xmax,Ymin,Ymax,Zmin,Zmax; + int getBrickId (double X, double Y, double Z) + { + int Ix = (int)((double)Nx * (X-Xmin) / (Xmax-Xmin)); + int Iy = (int)((double)Ny * (Y-Ymin) / (Ymax-Ymin)); + int Iz = (int)((double)Nz * (Z-Zmin) / (Zmax-Zmin)); + Ix = MIN(Ix,Nx-1); + Iy = MIN(Iy,Ny-1); + Iz = MIN(Iz,Nz-1); + if(Ix<0)Ix=0; + if(Iy<0)Iy=0; + if(Iz<0)Iz=0; + int index = Ix + Iy * Nx + Iz * Nx * Ny; + return index; + } + Brick<T> * getBrick (int index) + { + if(index <0 || index >= Nx*Ny*Nz)throw index; + Brick<T> *b = &(theTable[index]); + return b; + } + GetBox getBox; + InteriorCheck interiorCheck; + public: + mGeomSearch (double x1, + double x2, + double y1, + double y2, + double z1, + double z2, + GetBox g, + InteriorCheck i, + int nx = 10, + int ny = 10, + int nz = 10) : Nx(nx), Ny(ny), Nz(nz),getBox(g),interiorCheck(i) + { + Xmin = x1-TOL; Xmax = x2+TOL; + Ymin = y1-TOL; Ymax = y2+TOL; + Zmin = z1-TOL; Zmax = z2+TOL; + + theTable = new Brick<T> [Nx*Ny*Nz]; + } + + ~mGeomSearch () + { + delete [] theTable; + }; + Brick<T> *getBrick(double X, double Y, double Z) + { + return (getBrick(getBrickId(X,Y,Z))); + } + T * find (double X[3] , double U[3]); + bool AddObject ( T *); + bool RemoveObject ( T *); + }; + + template <class T,class GetBox, class InteriorCheck> + T * mGeomSearch<T,GetBox,InteriorCheck> :: find ( double X[3] , + double U[3] ) + { + Brick<T> *b = getBrick (X[0],X[1],X[2]); + if(!b) return 0; + + for(int i=0;i<b->size();i++) + { + T * obj = (*b)[i]; + if(interiorCheck(obj,X,U))return obj; + } + return 0; + } + + template <class T,class GetBox, class InteriorCheck> + bool mGeomSearch<T,GetBox, InteriorCheck> :: AddObject ( T * obj ) + { + int Ix1,Ix2,Iy1,Iy2,Iz1,Iz2; + int i,j,k,index; + Brick<T> *pBrick; + + /*Template Objects must overload getBox function*/ + double min[3]; + double max[3]; + getBox(obj,min,max); + + Ix1 = (int)( (double)Nx * (min[0] - Xmin) /( Xmax - Xmin )); + Ix2 = (int)( (double)Nx * (max[0] - Xmin) /( Xmax - Xmin )); + Iy1 = (int)( (double)Ny * (min[1] - Ymin) /( Ymax - Ymin )); + Iy2 = (int)( (double)Ny * (max[1] - Ymin) /( Ymax - Ymin )); + Iz1 = (int)( (double)Nz * (min[2] - Zmin) /( Zmax - Zmin )); + Iz2 = (int)( (double)Nz * (max[2] - Zmin) /( Zmax - Zmin )); + + + Ix1 = MAX(Ix1,0); + Ix2 = MIN(Ix2,Nx-1); + Iy1 = MAX(Iy1,0); + Iy2 = MIN(Iy2,Ny-1); + Iz1 = MAX(Iz1,0); + Iz2 = MIN(Iz2,Nz-1); + + for(i=Ix1;i<=Ix2;i++){ + for(j=Iy1;j<=Iy2;j++){ + for(k=Iz1;k<=Iz2;k++){ + index = i + j * Nx + k * Nx * Ny; + pBrick = getBrick(index); + pBrick->Objects.push_back(obj); + } + } + } + return true; + } + + template <class T,class GetBox, class InteriorCheck> + bool mGeomSearch<T,GetBox, InteriorCheck> :: RemoveObject ( T * obj ) + { + int Ix1,Ix2,Iy1,Iy2,Iz1,Iz2; + int i,j,k,index; + Brick<T> *pBrick; + + /*Template Objects must overload getBox function*/ + double min[3]; + double max[3]; + getBox(obj,min,max); + + Ix1 = (int)( (double)Nx * (min[0] - Xmin) /( Xmax - Xmin )); + Ix2 = (int)( (double)Nx * (max[0] - Xmin) /( Xmax - Xmin )); + Iy1 = (int)( (double)Ny * (min[1] - Ymin) /( Ymax - Ymin )); + Iy2 = (int)( (double)Ny * (max[1] - Ymin) /( Ymax - Ymin )); + Iz1 = (int)( (double)Nz * (min[2] - Zmin) /( Zmax - Zmin )); + Iz2 = (int)( (double)Nz * (max[2] - Zmin) /( Zmax - Zmin )); + + + Ix1 = MAX(Ix1,0); + Ix2 = MIN(Ix2,Nx-1); + Iy1 = MAX(Iy1,0); + Iy2 = MIN(Iy2,Ny-1); + Iz1 = MAX(Iz1,0); + Iz2 = MIN(Iz2,Nz-1); + + for(i=Ix1;i<=Ix2;i++){ + for(j=Iy1;j<=Iy2;j++){ + for(k=Iz1;k<=Iz2;k++){ + index = i + j * Nx + k * Nx * Ny; + pBrick = getBrick(index); + pBrick->Objects.erase ( std::remove (pBrick->Objects.begin(),pBrick->Objects.end(),obj) , + pBrick->Objects.end () ); + } + } + } + return true; + } + +} // end of namespace + +#undef TOL +#endif diff --git a/Plugin/Makefile b/Plugin/Makefile index 85a28768ed..7c9ea138de 100644 --- a/Plugin/Makefile +++ b/Plugin/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.20 2001-10-30 09:52:40 geuzaine Exp $ +# $Id: Makefile,v 1.21 2001-10-31 16:33:46 remacle Exp $ # # Makefile for "libAdapt.a" # @@ -18,7 +18,7 @@ OS_FLAGS = VERSION_FLAGS = RMFLAGS = -f -CFLAGS = $(OPT_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE) $(GUI_INCLUDE) +CFLAGS = -DMPICH_SKIP_MPICXX $(C_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE) $(GUI_INCLUDE) SRC = Plugin.cpp\ LevelsetPlugin.cpp\ diff --git a/Plugin/Plugin.cpp b/Plugin/Plugin.cpp index affcc20ed6..f995efd9f1 100644 --- a/Plugin/Plugin.cpp +++ b/Plugin/Plugin.cpp @@ -1,5 +1,6 @@ -// $Id: Plugin.cpp,v 1.21 2001-10-25 07:22:46 geuzaine Exp $ +// $Id: Plugin.cpp,v 1.22 2001-10-31 16:33:46 remacle Exp $ +#include <map> #ifndef _NODLL #include <dlfcn.h> #endif diff --git a/Triangle/Makefile b/Triangle/Makefile index f583d2cbbc..25d23f8800 100644 --- a/Triangle/Makefile +++ b/Triangle/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.3 2001-08-20 08:25:24 geuzaine Exp $ +# $Id: Makefile,v 1.4 2001-10-31 16:33:46 remacle Exp $ # # Makefile for "libTriangle.a" # @@ -24,7 +24,7 @@ .IGNORE: -CC = cc +CC = gcc AR = ar ruvs RM = rm RANLIB = ranlib diff --git a/benchmarks/3d/magnetron/magnetron1-full.geo b/benchmarks/3d/magnetron/magnetron1-full.geo index 1272164865..bb92a93228 100644 --- a/benchmarks/3d/magnetron/magnetron1-full.geo +++ b/benchmarks/3d/magnetron/magnetron1-full.geo @@ -1,5 +1,5 @@ mm = 0.001 ; // 1 milimetre = 0.001 metre -lc = 2.1*mm ; // unite de base min pour la taille caracteristique du maillage +lc = 0.75*mm ; // unite de base min pour la taille caracteristique du maillage lcpba2 = 1.5*lc ; // lc dessous plaque base lcpba1 = lc ; // lc dessus plaque base -- GitLab