1 // ********** DO NOT REMOVE THIS BANNER ********** 2 // ORIG-DATE: Jan 2008 3 // -*- Mode : c++ -*- 4 // 5 // SUMMARY : Generic Tree header (binairy, Quad, Oct) 6 // USAGE : LGPL 7 // ORG : LJLL Universite Pierre et Marie Curi, Paris, FRANCE 8 // AUTHOR : Frederic Hecht 9 // E-MAIL : frederic.hecht@ann.jussieu.fr 10 // 11 12 /* 13 14 This file is part of Freefem++ 15 16 Freefem++ is free software; you can redistribute it and/or modify 17 it under the terms of the GNU Lesser General Public License as published by 18 the Free Software Foundation; either version 2.1 of the License, or 19 (at your option) any later version. 20 21 Freefem++ is distributed in the hope that it will be useful, 22 but WITHOUT ANY WARRANTY; without even the implied warranty of 23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 24 GNU Lesser General Public License for more details. 25 26 You should have received a copy of the GNU Lesser General Public License 27 along with Freefem++; if not, write to the Free Software 28 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 29 30 Thank to the ARN () FF2A3 grant 31 ref:ANR-07-CIS7-002-01 32 */ 33 34 namespace Fem2D { 35 #include "R3.hpp" 36 } 37 namespace EF23 { 38 using namespace Fem2D; 39 using Fem2D::R; 40 41 static const int MaxDeep = 30; 42 typedef int IntQuad; 43 typedef long long Int8; 44 static const IntQuad MaxISize = ( 1L << MaxDeep); 45 static const IntQuad MaxISize1 = MaxISize-1; 46 class Z1 { public: INTER_SEG1d(int a,int b,int x,int y)47 static bool INTER_SEG1d(int a,int b,int x,int y) { return (((y) > (a)) && ((x) <(b)));} 48 int x; Z1()49 Z1():x(0){} Z1(R1 P)50 Z1(R1 P) : x((int)P.x) {} Z1(int i)51 Z1( int i) : x(i){} Z1(const Z1 & pp,int k,int l)52 Z1(const Z1 &pp,int k,int l): x(pp.x+(( k&1) ? l : 0)) {} Add(int k,int l)53 void Add(int k,int l) { x+= (( k&1) ? l : 0) ;} Z1(const Z1 & A,const Z1 & B)54 Z1(const Z1 &A,const Z1 &B) : x(B.x-A.x){} 55 Case(int l) const56 int Case(int l) const { return ( x & l)? 1 : 0 ;} norm() const57 int norm() const { return abs(x);} norm2() const58 Int8 norm2() const { return (Int8) x* (Int8) x;} Bound()59 void Bound() { x = max(min(x,MaxISize1),0);} 60 less(Z1 h) const61 bool less(Z1 h) const { return abs(x) <h.x ;} interseg(Z1 pp,int hb,int h) const62 bool interseg(Z1 pp,int hb,int h) const { 63 return INTER_SEG1d(x,x+hb,pp.x-h,pp.x+h) ;} interseg(const Z1 & pp,int hb,const Z1 & h) const64 bool interseg(const Z1 &pp,int hb,const Z1 &h) const { 65 return INTER_SEG1d(x,x+hb,pp.x-h.x,pp.x+h.x) ;} operator R1() const66 operator R1 () const { return R1(x);} 67 }; 68 69 70 71 class Z2 { public: INTER_SEG1d(int a,int b,int x,int y)72 static bool INTER_SEG1d(int a,int b,int x,int y) { return (((y) > (a)) && ((x) <(b)));} 73 int x,y; Z2()74 Z2():x(0),y(0) {} Z2(R2 P)75 Z2(R2 P) : x((int)P.x),y((int)P.y) {} Z2(int i)76 Z2( int i) : x(i),y(i){} 77 //Z2( int i,int j) : x(i),y(j) {} Z2(const Z2 & pp,int k,int l)78 Z2(const Z2 &pp,int k,int l): x(pp.x+(( k&1) ? l : 0)),y(pp.y+(( k&2) ? l : 0)) {} Add(int k,int l)79 void Add(int k,int l) { x+= (( k&1) ? l : 0) ; y+= (( k&2) ? l : 0);} Z2(const Z2 & A,const Z2 & B)80 Z2(const Z2 &A,const Z2 &B) : x(B.x-A.x),y(B.y-A.y) {} 81 Case(int l) const82 int Case(int l) const { return ( ( y & l) ? (( x & l) ? 3 : 2 ) : ( ( x & l)? 1 : 0 ) ) ;} norm() const83 int norm() const { return Max(abs(x),abs(y));} norm2() const84 Int8 norm2() const { return (Int8) x*(Int8) x + (Int8) y*(Int8)y;} Bound()85 void Bound() { x = max(min(x,MaxISize1),0); 86 y = max(min(y,MaxISize1),0);} 87 less(Z2 h) const88 bool less(Z2 h) const { return abs(x) <h.x && abs(y) <h.y;} interseg(Z2 pp,int hb,int h) const89 bool interseg(Z2 pp,int hb,int h) const { 90 return INTER_SEG1d(x,x+hb,pp.x-h,pp.x+h) && INTER_SEG1d(y,y+hb,pp.y-h,pp.y+h);} interseg(const Z2 & pp,int hb,const Z2 & h) const91 bool interseg(const Z2 &pp,int hb,const Z2 &h) const { 92 return INTER_SEG1d(x,x+hb,pp.x-h.x,pp.x+h.x) && INTER_SEG1d(y,y+hb,pp.y-h.y,pp.y+h.y);} operator R2() const93 operator R2 () const { return R2(x,y);} 94 }; 95 96 class Z3 { public: INTER_SEG1d(int a,int b,int x,int y)97 static bool INTER_SEG1d(int a,int b,int x,int y) { return (((y) > (a)) && ((x) <(b)));} 98 int x,y,z; Z3()99 Z3():x(0),y(0),z(0) {} 100 Z3(R3 P)101 Z3(R3 P) : x((int)P.x),y((int)P.y),z((int) P.z) {} Z3(int i)102 Z3( int i) : x(i),y(i),z(i){} 103 Z3(const Z3 & pp,int k,int l)104 Z3(const Z3 &pp,int k,int l): x(pp.x+(( k&1) ? l : 0)),y(pp.y+(( k&2) ? l : 0)),z(pp.z+(( k&4) ? l : 0)) {} Add(int k,int l)105 void Add(int k,int l) { x+= (( k&1) ? l : 0) ; y+= (( k&2) ? l : 0); z+= (( k&4) ? l : 0);} Z3(const Z3 & A,const Z3 & B)106 Z3(const Z3 &A,const Z3 &B) : x(B.x-A.x),y(B.y-A.y),z(B.z-A.z) {} Bound()107 void Bound() { x = max(min(x,MaxISize1),0); 108 y = max(min(y,MaxISize1),0); 109 z = max(min(z,MaxISize1),0);} 110 Case(int l) const111 int Case(int l) const {// cout << " case = "<< int((x&l)!=0)+(int((y&l)!=0)<<1) + (int((z&l)!=0)<<2); 112 return int( (x&l)!=0) + ( int((y&l)!=0)<<1 ) + ( int( (z&l)!=0) <<2 ) ;} norm() const113 int norm() const { return Max(abs(x),abs(y),abs(z));} norm2() const114 Int8 norm2() const { return (Int8) x*(Int8) x + (Int8) y*(Int8)y + (Int8) z*(Int8)z; } less(Z3 h) const115 bool less(Z3 h) const { return abs(x) <h.x && abs(y) <h.y && abs(z) < h.z ;} interseg(Z3 pp,int hb,int h) const116 bool interseg(Z3 pp,int hb,int h) const { 117 return INTER_SEG1d(x,x+hb,pp.x-h,pp.x+h) && INTER_SEG1d(y,y+hb,pp.y-h,pp.y+h) && INTER_SEG1d(z,z+hb,pp.z-h,pp.z+h) ; 118 } interseg(const Z3 & pp,int hb,const Z3 & h) const119 bool interseg(const Z3 &pp,int hb,const Z3 &h) const { 120 return INTER_SEG1d(x,x+hb,pp.x-h.x,pp.x+h.x) && INTER_SEG1d(y,y+hb,pp.y-h.y,pp.y+h.y) && INTER_SEG1d(z,z+hb,pp.z-h.z,pp.z+h.z); 121 } operator R3() const122 operator R3 () const { return R3(x,y,z);} 123 124 }; 125 operator <<(ostream & f,const Z3 & P)126 inline ostream& operator <<(ostream& f, const Z3 & P ) { f << P.x << ' ' << P.y << ' ' << P.z ; return f; } operator <<(ostream & f,const Z2 & P)127 inline ostream& operator <<(ostream& f, const Z2 & P ) { f << P.x << ' ' << P.y ; return f; } operator <<(ostream & f,const Z1 & P)128 inline ostream& operator <<(ostream& f, const Z1 & P ) { f << P.x ; return f; } 129 130 template<class Rd> struct Traits_Zd { typedef void Zd;}; 131 template<> struct Traits_Zd<R1> { typedef Z1 Zd;}; 132 template<> struct Traits_Zd<R2> { typedef Z2 Zd;}; 133 template<> struct Traits_Zd<R3> { typedef Z3 Zd;}; 134 135 template<class Vertex> 136 class GTree { 137 typedef typename Vertex::Rd Rd; 138 typedef typename Traits_Zd<Rd>::Zd Zd; 139 140 141 public: 142 143 static const int d =Rd::d; 144 static const int N = 1 << d; // N=2^(d-1) 145 146 147 class QuadTreeBox { 148 public: 149 150 int n; // if n < 4 => Vertex else => QuadTreeBox; 151 union { 152 QuadTreeBox *b[N]; 153 Vertex * v[N]; 154 }; 155 // void init() { for(int i=0;i<N;++i) b[i]=0;} 156 157 }; // end class QuadTreeBox ///////////////// 158 159 class StorageQuadTreeBox { 160 public: 161 QuadTreeBox *b,*bc,*be; 162 int len; 163 StorageQuadTreeBox *n; // next StorageQuadTreeBox 164 StorageQuadTreeBox(int ,StorageQuadTreeBox * =0); 165 ~StorageQuadTreeBox(); SizeOf() const166 int SizeOf() const { 167 return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0); 168 } 169 }; // end class StorageQuadTreeBox 170 171 StorageQuadTreeBox * sb; 172 173 174 int lenStorageQuadTreeBox; 175 176 public: 177 QuadTreeBox * root; 178 // Mesh *th; 179 180 int NbQuadTreeBoxSearch,NbVerticesSearch; 181 int NbQuadTreeBox,NbVertices; 182 183 Rd cMin,cMax; // box of QuadTree 184 R coef; // 185 186 RdtoZd(const Rd & P) const187 Zd RdtoZd(const Rd &P) const {return Zd((Minc(Maxc(P,cMin),cMax)-cMin)*coef);} VtoZd(const Vertex * v) const188 Zd VtoZd(const Vertex * v) const {return RdtoZd( (const Rd&) *v);} VtoZd(const Vertex & v) const189 Zd VtoZd(const Vertex & v) const {return RdtoZd( (const Rd&) v);} 190 ZdtoRd(const Zd & I) const191 Rd ZdtoRd(const Zd &I) const { return ( (Rd) I )/coef+cMin;} 192 NearestVertex(const Rd & P,bool trueNearest=false)193 Vertex * NearestVertex(const Rd & P,bool trueNearest=false) { 194 return NearestVertex(RdtoZd(P),trueNearest);} //XtoI(P.x),YtoJ(P.y));} 195 Vertex * NearestVertexWithNormal(const Rd & P); 196 Vertex * NearestVertex(Zd i2,bool trueNearest=false); 197 198 Vertex * ToClose(const Rd & ,R ,Zd, bool nearest=false ); ToClose(const Rd & P,R delta,bool nearest=false)199 Vertex * ToClose(const Rd & P,R delta,bool nearest=false){ 200 int hx = (int) (coef*delta); 201 hx = hx>0 ? hx:1; // bof bof .... 202 // ffassert(hx>0);// bug if too small dec. 2019 FH. 203 //if(verbosity > 5 ) cout << "hx=" << hx << " coef=" << coef << endl; 204 return ToClose(P,delta,Zd(hx),nearest);} SizeOf() const205 int SizeOf() const {return sizeof(GTree)+sb->SizeOf();} 206 207 void Add( Vertex & w); 208 NewQuadTreeBox()209 QuadTreeBox* NewQuadTreeBox() 210 { 211 ///cout << "NewQuadTreeBox " << sb << " " << sb->bc << " " 212 //<< sb->be << " " <<lenStorageQuadTreeBox <<endl; 213 if(! (sb->bc<sb->be)) 214 sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); 215 216 assert(sb && (sb->bc->n == 0)); 217 NbQuadTreeBox++; 218 return sb->bc++; 219 } 220 ~GTree(); 221 GTree(Vertex * v,Rd Pmin,Rd Pmax,int nbv); 222 GTree(); 223 template<class V> 224 friend ostream& operator <<(ostream& f, const GTree<V> & qt); 225 // add FH mars 2020 , for new search .. FH... ListNearestVertex(Vertex ** lnv,int nvn,double delta,Rd P)226 int ListNearestVertex(Vertex **lnv,int nvn,double delta,Rd P) 227 { int hx = (int) (coef*delta); 228 hx = hx>0 ? hx:1; // bof bof .... 229 return ListNearestVertex(lnv,nvn,hx,RdtoZd(P));} TrueNearestVertex(Zd i2)230 Vertex * TrueNearestVertex(Zd i2) {return NearestVertex(i2,true);} 231 int ListNearestVertex(Vertex **lnv,int nlvnx,int dh,Zd xyi); 232 233 TrueNearestVertex(const Rd & P)234 Vertex * TrueNearestVertex(const Rd & P) { 235 return TrueNearestVertex(RdtoZd(P));} 236 237 }; 238 239 template<class Mesh> 240 const typename Mesh::Element * Find(const Mesh & Th, 241 GTree<typename Mesh::Vertex> *quadtree, 242 typename Mesh::Rd P, 243 typename Mesh::RdHat & Phat, 244 bool & outside, 245 const typename Mesh::Element * tstart); 246 247 248 } // name space 249 250