/[escript]/branches/refine/buckley/src/OctCell1.cc
ViewVC logotype

Contents of /branches/refine/buckley/src/OctCell1.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3668 - (show annotations)
Wed Nov 16 01:49:46 2011 UTC (7 years, 10 months ago) by jfenwick
File size: 9948 byte(s)
Stage 1 rename
1 #include <iostream>
2 #include "OctCell.h"
3
4 using namespace buckley;
5
6 using namespace std;
7
8 namespace
9 {
10 const double eps=1e-1;
11
12 }
13
14 OctCell::OctCell(double cx, double cy, double cz, double wx, double wy, double wz, OctCell* progenitor)
15 {
16 leaf=true;
17 centre[0]=cx;
18 centre[1]=cy;
19 centre[2]=cz;
20 sides[0]=wx;
21 sides[1]=wy;
22 sides[2]=wz;
23 parent=progenitor;
24 if (parent)
25 {
26 depth=parent->depth+1;
27 }
28 else
29 {
30 depth=0;
31 }
32 id=0;
33 }
34
35 OctCell::~OctCell()
36 {
37 if (!leaf)
38 {
39 for (unsigned int i=0;i<8;++i)
40 {
41 delete kids[i];
42 }
43 }
44 }
45
46 // If we are allocating from a pool we need to keep a record of it somewhere
47 // maybe keep track (static?) of holding object with operators
48 // or define OctCell::new ?
49 void OctCell::split()
50 {
51 if (leaf)
52 {
53 const double dx=sides[0]/2;
54 const double dy=sides[1]/2;
55 const double dz=sides[2]/2;
56 const double cx=centre[0];
57 const double cy=centre[1];
58 const double cz=centre[2];
59 kids[0]=new OctCell(cx-dx/2, cy-dy/2, cz-dz/2, dx, dy, dz, this);
60 kids[1]=new OctCell(cx+dx/2, cy-dy/2, cz-dz/2, dx, dy, dz, this);
61 kids[2]=new OctCell(cx+dx/2, cy+dy/2, cz-dz/2, dx, dy, dz, this);
62 kids[3]=new OctCell(cx-dx/2, cy+dy/2, cz-dz/2, dx, dy, dz, this);
63 kids[4]=new OctCell(cx-dx/2, cy-dy/2, cz+dz/2, dx, dy, dz, this);
64 kids[5]=new OctCell(cx+dx/2, cy-dy/2, cz+dz/2, dx, dy, dz, this);
65 kids[6]=new OctCell(cx+dx/2, cy+dy/2, cz+dz/2, dx, dy, dz, this);
66 kids[7]=new OctCell(cx-dx/2, cy+dy/2, cz+dz/2, dx, dy, dz, this);
67 leaf=false;
68 }
69 }
70
71 // could use the replacement delete operator to notify other parts of the system that a value has been removed.
72 void OctCell::merge()
73 {
74 if (!leaf)
75 {
76 for (unsigned int i=0;i<8;i++)
77 {
78 delete kids[i];
79 }
80 leaf=true;
81 }
82 }
83
84 void OctCell::allSplit(unsigned int depth)
85 {
86 if (depth>0)
87 {
88 if (leaf)
89 {
90 split();
91 }
92 if (depth==1)
93 {
94 return;
95 }
96 for (unsigned int i=0;i<8;i++)
97 {
98 kids[i]->split();
99 kids[i]->allSplit(depth-1);
100 }
101 }
102 }
103
104
105 // After a cell has been split, this method can be called to ensure that the tree is still "not too unbalanced"
106 void OctCell::outwardRefine(unsigned desireddepth)
107 {
108 double minside=(sides[0]>sides[1])?sides[0]:sides[1];
109 minside=(minside>sides[2])?sides[2]:minside;
110 minside/=5; // this will work provided that all other divisions have been safe
111
112 // cerr << "Outward buckley: " << centre[0] << " " << centre[1] << " " << centre[2] << "\n";
113 // cerr << "-\n";
114 upSplitPoint(centre[0]-sides[0]/2-minside, centre[1], centre[2], desireddepth);
115 // cerr << "-\n";
116 upSplitPoint(centre[0]+sides[0]/2+minside, centre[1], centre[2], desireddepth);
117 // cerr << "-\n";
118 upSplitPoint(centre[0], centre[1]-sides[1]/2-minside, centre[2], desireddepth);
119 // cerr << "-\n";
120 upSplitPoint(centre[0], centre[1]+sides[1]/2+minside, centre[2], desireddepth);
121 // cerr << "-\n";
122 upSplitPoint(centre[0], centre[1], centre[2]-sides[2]/2-minside, desireddepth);
123 // cerr << "-\n";
124 upSplitPoint(centre[0], centre[1], centre[2]+sides[2]/2+minside, desireddepth);
125 // cerr << "----\n";
126 }
127
128 // Works up the tree until it finds
129 void OctCell::upSplitPoint(double x, double y, double z, unsigned d)
130 {
131 if ((x<centre[0]-sides[0]/2) || (x>centre[0]+sides[0]/2) || (y<centre[1]-sides[1]/2) ||
132 (y>centre[1]+sides[1]/2) || (z<centre[2]-sides[2]/2) || (z>centre[2]+sides[2]/2))
133 {
134 // It is outside this box so go up
135 if (parent!=0)
136 {
137 parent->upSplitPoint(x, y, z, d);
138 }
139 return;
140 }
141 else
142 {
143 // cerr << "Up to " << centre[0] << " " << centre[1] << " " << centre[2] << "\n";
144
145 // it's in this cell somewhere so start moving down again
146 splitPoint(x, y, z, d);
147 }
148 }
149
150 // Works up the tree until it finds
151 void OctCell::upCollPoint(double x, double y, double z, unsigned d)
152 {
153 if ((x<centre[0]-sides[0]/2) || (x>centre[0]+sides[0]/2) || (y<centre[1]-sides[1]/2) ||
154 (y>centre[1]+sides[1]/2) || (z<centre[2]-sides[2]/2) || (z>centre[2]+sides[2]/2))
155 {
156 // It is outside this box so go up
157 if (parent!=0)
158 {
159 parent->upCollPoint(x, y, z, d);
160 }
161 return;
162 }
163 else
164 {
165 // cerr << "Up to " << centre[0] << " " << centre[1] << " " << centre[2] << "\n";
166
167 // it's in this cell somewhere so start moving down again
168 collapsePoint(x, y, z, d);
169 }
170 }
171
172
173 // return the leaf which contains this point.
174 // Assumes that the point lies somewhere in this cell
175 OctCell* OctCell::findLeaf(double x, double y, double z)
176 {
177 if (leaf)
178 {
179 return this;
180 }
181 int pz=(z>centre[2])?4:0;
182 int py=(y>centre[1]);
183 int px=(x>centre[0]);
184 int child=pz+2*py+(px^py);
185 return kids[child]->findLeaf(x,y,z);
186 }
187
188
189 // divide the leaf containing the point (if required)
190 // This needs to be more efficient but will do for now
191 void OctCell::splitPoint(double x, double y, double z, unsigned d)
192 {
193 // find cell
194 // if it doesn't need refining, then bail
195 // buckley neighbours to be at least d-1
196 // buckley this cell
197 OctCell* start=this;
198 do
199 {
200 OctCell* l=start->findLeaf(x, y, z);
201 if (l->depth>=d)
202 {
203 return;
204 }
205 l->outwardRefine(l->depth);
206 l->split();
207 start=l;
208 } while (start->depth+1 <d);
209 }
210
211
212 // collapses the children of this node
213 void OctCell::collapse()
214 {
215 if (!leaf)
216 {
217 // collapse any grandkids
218 for (unsigned i=0;i<8;++i)
219 {
220 kids[i]->collapse(); // so we know we are only dealing with leaves
221 }
222 for (unsigned i=0;i<8;++i)
223 {
224 delete kids[i];
225 }
226 leaf=true;
227 }
228 }
229
230 void OctCell::collapseAll(unsigned int desdepth)
231 {
232 if (leaf)
233 {
234 return;
235 }
236 for (unsigned i=0;i<8;++i)
237 {
238 kids[i]->collapseAll(desdepth);
239 }
240 if (depth>=desdepth)
241 {
242 collapse();
243 }
244 }
245
246
247 // The leaf which contains (x,y,z) should be at depth at most d
248 void OctCell::collapsePoint(double x, double y, double z, unsigned d)
249 {
250 // find cell
251 // if it doesn't need refining, then bail
252 // buckley neighbours to be at least d-1
253 // buckley this cell
254 OctCell* start=this;
255 do
256 {
257 OctCell* l=start->findLeaf(x, y, z);
258 if (l->depth<=d)
259 {
260 return;
261 }
262 l=l->parent;
263 l->outwardCollapse(l->depth+1); // since we have gone up a level
264 l->collapse();
265 start=l;
266 } while (start->depth+1 >d);
267
268 }
269
270
271 // This is horribly inefficient!!!
272 // After a cell has been, thisdegrem method can be called to ensure that the tree is still "not too unbalanced"
273 void OctCell::outwardCollapse(unsigned desireddepth)
274 {
275 double minside=(sides[0]>sides[1])?sides[0]:sides[1];
276 minside=(minside>sides[2])?sides[2]:minside;
277 minside/=5; // this will work provided that all other divisions have been safe
278
279 // cerr << "Outward collapse: " << centre[0] << " " << centre[1] << " " << centre[2] << "\n";
280 // cerr << "-" << (centre[0]-sides[0]/2-minside) << ' ' << centre[1] << ' ' << centre[2]<<"\n";
281 upCollPoint(centre[0]-sides[0]/2-minside, centre[1]-minside, centre[2]-minside, desireddepth);
282
283 upCollPoint(centre[0]-sides[0]/2-minside, centre[1]-minside, centre[2]+minside, desireddepth);
284
285 upCollPoint(centre[0]-sides[0]/2-minside, centre[1]+minside, centre[2]-minside, desireddepth);
286
287 upCollPoint(centre[0]-sides[0]/2-minside, centre[1]+minside, centre[2]+minside, desireddepth);
288
289
290 // cerr << "-\n";
291 upCollPoint(centre[0]+sides[0]/2+minside, centre[1]-minside, centre[2]-minside, desireddepth);
292
293 upCollPoint(centre[0]+sides[0]/2+minside, centre[1]-minside, centre[2]+minside, desireddepth);
294
295 upCollPoint(centre[0]+sides[0]/2+minside, centre[1]+minside, centre[2]-minside, desireddepth);
296
297 upCollPoint(centre[0]+sides[0]/2+minside, centre[1]+minside, centre[2]+minside, desireddepth);
298
299
300 // cerr << "-\n";
301 upCollPoint(centre[0]-minside, centre[1]-sides[1]/2-minside, centre[2]-minside, desireddepth);
302
303 upCollPoint(centre[0]-minside, centre[1]-sides[1]/2-minside, centre[2]+minside, desireddepth);
304
305 upCollPoint(centre[0]+minside, centre[1]-sides[1]/2-minside, centre[2]-minside, desireddepth);
306
307 upCollPoint(centre[0]+minside, centre[1]-sides[1]/2-minside, centre[2]+minside, desireddepth);
308
309
310
311
312
313 // cerr << "-\n";
314 upCollPoint(centre[0]-minside, centre[1]+sides[1]/2+minside, centre[2]-minside, desireddepth);
315
316 upCollPoint(centre[0]-minside, centre[1]+sides[1]/2+minside, centre[2]+minside, desireddepth);
317
318 upCollPoint(centre[0]+minside, centre[1]+sides[1]/2+minside, centre[2]-minside, desireddepth);
319
320 upCollPoint(centre[0]+minside, centre[1]+sides[1]/2+minside, centre[2]+minside, desireddepth);
321
322
323 // cerr << "-\n";
324 upCollPoint(centre[0]-minside, centre[1]-minside, centre[2]-sides[2]/2-minside, desireddepth);
325
326 upCollPoint(centre[0]-minside, centre[1]+minside, centre[2]-sides[2]/2-minside, desireddepth);
327
328 upCollPoint(centre[0]+minside, centre[1]-minside, centre[2]-sides[2]/2-minside, desireddepth);
329
330 upCollPoint(centre[0]+minside, centre[1]+minside, centre[2]-sides[2]/2-minside, desireddepth);
331
332 // cerr << "-\n";
333 upCollPoint(centre[0]-minside, centre[1]-minside, centre[2]+sides[2]/2+minside, desireddepth);
334
335 upCollPoint(centre[0]-minside, centre[1]+minside, centre[2]+sides[2]/2+minside, desireddepth);
336
337 upCollPoint(centre[0]+minside, centre[1]-minside, centre[2]+sides[2]/2+minside, desireddepth);
338
339 upCollPoint(centre[0]+minside, centre[1]+minside, centre[2]+sides[2]/2+minside, desireddepth);
340
341 // cerr << "----\n";
342 }
343
344 void OctCell::doLeafWalk(cellfunct c, void* v)
345 {
346 if (leaf)
347 {
348 c(*this, v);
349 }
350 else
351 {
352 for (unsigned int i=0;i<8;++i)
353 {
354 kids[i]->doLeafWalk(c,v);
355 }
356 }
357 }
358
359
360
361

  ViewVC Help
Powered by ViewVC 1.1.26