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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3680 - (show annotations)
Fri Nov 18 04:48:53 2011 UTC (7 years, 10 months ago) by jfenwick
File size: 11995 byte(s)
It has memory errors during a getX.

1 #include <iostream>
2 #include "OctCell.h"
3 #include "LeafInfo.h"
4
5 using namespace buckley;
6
7 using namespace std;
8
9 namespace
10 {
11 const double eps=1e-1;
12 }
13
14
15 #include "FaceConsts.h"
16
17 OctCell::OctCell(double cx, double cy, double cz, double wx, double wy, double wz, OctCell* progenitor)
18 {
19 leaf=true;
20 centre[0]=cx;
21 centre[1]=cy;
22 centre[2]=cz;
23 sides[0]=wx;
24 sides[1]=wy;
25 sides[2]=wz;
26 parent=progenitor;
27 if (parent)
28 {
29 depth=parent->depth+1;
30 }
31 else
32 {
33 depth=0;
34 }
35 id=0;
36 leafinfo=new LeafInfo(this);
37 }
38
39 OctCell::~OctCell()
40 {
41 if (!leaf)
42 {
43 for (unsigned int i=0;i<8;++i)
44 {
45 delete kids[i];
46 }
47 }
48 else
49 {
50 delete leafinfo;
51 }
52 }
53
54
55
56
57
58 // If we are allocating from a pool we need to keep a record of it somewhere
59 // maybe keep track (static?) of holding object with operators
60 // or define OctCell::new ?
61 void OctCell::split(size_t* leafc)
62 {
63 if (leaf)
64 {
65 for (int f=0;f<6;f++)
66 {
67 if (leafinfo->next[f] && (depth>leafinfo->next[f]->owner->depth))
68 {
69 leafinfo->next[f]->owner->split(leafc);
70 }
71
72 }
73 const double dx=sides[0]/2;
74 const double dy=sides[1]/2;
75 const double dz=sides[2]/2;
76 const double cx=centre[0];
77 const double cy=centre[1];
78 const double cz=centre[2];
79 kids[0]=new OctCell(cx-dx/2, cy-dy/2, cz-dz/2, dx, dy, dz, this);
80 kids[1]=new OctCell(cx+dx/2, cy-dy/2, cz-dz/2, dx, dy, dz, this);
81 kids[2]=new OctCell(cx+dx/2, cy+dy/2, cz-dz/2, dx, dy, dz, this);
82 kids[3]=new OctCell(cx-dx/2, cy+dy/2, cz-dz/2, dx, dy, dz, this);
83 kids[4]=new OctCell(cx-dx/2, cy-dy/2, cz+dz/2, dx, dy, dz, this);
84 kids[5]=new OctCell(cx+dx/2, cy-dy/2, cz+dz/2, dx, dy, dz, this);
85 kids[6]=new OctCell(cx+dx/2, cy+dy/2, cz+dz/2, dx, dy, dz, this);
86 kids[7]=new OctCell(cx-dx/2, cy+dy/2, cz+dz/2, dx, dy, dz, this);
87 leaf=false;
88 leafinfo->split(kids);
89 delete leafinfo;
90 leafinfo=0;
91 (*leafc)+=7;
92 }
93 }
94
95 void OctCell::allSplit(unsigned int depth, size_t* leafc)
96 {
97 if (depth>0)
98 {
99 if (leaf)
100 {
101 split(leafc);
102 }
103 if (depth==1)
104 {
105 return;
106 }
107 for (unsigned int i=0;i<8;i++)
108 {
109 kids[i]->split(leafc);
110 kids[i]->allSplit(depth-1, leafc);
111 }
112 }
113 }
114
115 // return the leaf which contains this point.
116 // Assumes that the point lies somewhere in this cell
117 OctCell* OctCell::findLeaf(double x, double y, double z)
118 {
119 if (leaf)
120 {
121 return this;
122 }
123 int pz=(z>centre[2])?4:0;
124 int py=(y>centre[1]);
125 int px=(x>centre[0]);
126 int child=pz+2*py+(px^py);
127 return kids[child]->findLeaf(x,y,z);
128 }
129
130
131 // divide the leaf containing the point (if required)
132 // This needs to be more efficient but will do for now
133 void OctCell::splitPoint(double x, double y, double z, unsigned d, size_t* leafc)
134 {
135 // find cell
136 // if it doesn't need refining, then bail
137 // buckley neighbours to be at least d-1
138 // buckley this cell
139 OctCell* start=this;
140 do
141 {
142 OctCell* l=start->findLeaf(x, y, z);
143 if (l->depth>=d)
144 {
145 return;
146 }
147 l->split(leafc);
148 start=l;
149 } while (start->depth+1 <d);
150 }
151
152 // collapse the children of this cell into it
153 void OctCell::collapse(size_t* leafc)
154 {
155 if (!leaf)
156 {
157 // collapse any grandkids
158 for (unsigned i=0;i<8;++i)
159 {
160 if (! kids[i]->leaf)
161 {
162 kids[i]->collapse(leafc); // so we know we are only dealing with leaves
163 }
164 }
165 // now we need to ensure that collapsing this would not break any neighbours
166 /* // we do this by looking at faces from opposite corners
167
168 for (int j=1;j<7;j+=2)
169 {
170
171 if (kids[0]->leafinfo->next[j])
172 {
173 if (depth+1 <= kids[0]->leafinfo->next[j]->owner->depth)
174 {
175 if (depth+1<kids[0]->leafinfo->next[j]->owner->depth)
176 {
177 kids[0]->leafinfo->next[j]->owner->parent->collapse();
178 }
179 else // entering else means that there is at least one leaf on the face which is
180 { // at the correct level but that doesn't mean they all are
181 for (unsigned s=0;s<4;++s)
182 {
183 if (! kids[0]->leafinfo->next[j]->owner->parent->kids[faces[oppdir(j)][s]]->leaf)
184 {
185 kids[0]->leafinfo->next[j]->owner->parent->collapse();
186 }
187 }
188 }
189 }
190 }
191
192 }
193
194 for (int j=0;j<6;j+=2)
195 {
196
197 if (kids[6]->leafinfo->next[j])
198 {
199 if (depth+1 <= kids[6]->leafinfo->next[j]->owner->depth)
200 {
201 if (depth+1<kids[6]->leafinfo->next[j]->owner->depth)
202 {
203 kids[6]->leafinfo->next[j]->owner->parent->collapse();
204 }
205 else // entering else means that there is at least one leaf on the face which is
206 { // at the correct level but that doesn't mean they all are
207 for (unsigned s=0;s<4;++s)
208 {
209 if (! kids[6]->leafinfo->next[j]->owner->parent->kids[faces[oppdir(j)][s]]->leaf)
210 {
211 kids[6]->leafinfo->next[j]->owner->parent->collapse();
212 }
213 }
214 }
215 }
216 }
217
218 }*/
219
220 /*
221 if (kids[0]->leafinfo->next[1] && (depth+1<kids[0]->leafinfo->next[1]->owner->depth))
222 {
223 kids[0]->leafinfo->next[1]->owner->parent->collapse();
224 }
225 if (kids[0]->leafinfo->next[3] && (depth+1<kids[0]->leafinfo->next[3]->owner->depth))
226 {
227 kids[0]->leafinfo->next[3]->owner->parent->collapse();
228 }
229 if (kids[0]->leafinfo->next[5] && (depth+1<kids[0]->leafinfo->next[5]->owner->depth))
230 {
231 kids[0]->leafinfo->next[5]->owner->parent->collapse();
232 }*/
233 /*
234 if (kids[6]->leafinfo->next[0] && (depth+1<kids[6]->leafinfo->next[0]->owner->depth))
235 {
236 kids[6]->leafinfo->next[0]->owner->parent->collapse();
237 }
238 if (kids[6]->leafinfo->next[2] && (depth+1<kids[6]->leafinfo->next[2]->owner->depth))
239 {
240 kids[6]->leafinfo->next[2]->owner->parent->collapse();
241 }
242 if (kids[6]->leafinfo->next[4] && (depth+1<kids[6]->leafinfo->next[4]->owner->depth))
243 {
244 kids[6]->leafinfo->next[4]->owner->parent->collapse();
245 }*/
246 // }
247
248
249 unsigned newdepth=depth;
250 unsigned leafdepth=depth+1;
251
252 for (int k=0;k<6;++k)
253 {
254 int kid=faces[k][0];
255 LeafInfo* neigh=kids[kid]->leafinfo->next[k];
256
257 if (!neigh || (neigh->owner->depth < leafdepth)) // no neigbour or leaf is already shallower than our leaf
258 {
259 continue;
260 }
261 // now we need to check to see if there are any subdivided cells on this face
262 OctCell* nparent=neigh->owner->parent;
263 if (nparent->depth>newdepth) // this is in case we ended up pointing at a divided cell on this face
264 {
265 nparent=nparent->parent;
266 }
267 for (unsigned i=0;i<4;++i)
268 {
269 if (! nparent->kids[faces[oppdir(k)][i]]->leaf) // if there are any non-leaves on the neighbour face
270 {
271 nparent->kids[faces[oppdir(k)][i]]->collapse(leafc);
272 }
273 }
274 }
275 kids[0]->leafinfo->merge(); // pick a child and collapse it
276 leaf=true;
277 for (unsigned i=0;i<8;++i)
278 {
279 delete kids[i];
280 }
281 leaf=true;
282 (*leafc)-=7;
283 }
284 }
285
286 void OctCell::collapseAll(unsigned int desdepth, size_t* leafc)
287 {
288 if (leaf)
289 {
290 return;
291 }
292 for (unsigned i=0;i<8;++i)
293 {
294 kids[i]->collapseAll(desdepth, leafc);
295 }
296 if (depth>=desdepth)
297 {
298 collapse(leafc);
299 }
300 }
301
302
303 // The leaf which contains (x,y,z) should be at depth at most d
304 void OctCell::collapsePoint(double x, double y, double z, unsigned d, size_t* leafc)
305 {
306 // find cell
307 // if it doesn't need refining, then bail
308 // buckley neighbours to be at least d-1
309 // buckley this cell
310 OctCell* start=this;
311 do
312 {
313 OctCell* l=start->findLeaf(x, y, z);
314 if (l->depth<=d)
315 {
316 return;
317 }
318 l=l->parent;
319 l->collapse(leafc);
320 start=l;
321 } while (start->depth+1 >d);
322
323 }
324
325
326 void OctCell::doLeafWalk_const(const_cellfn c, void* v) const
327 {
328 if (leaf)
329 {
330 c(*this, v);
331 }
332 else
333 {
334 for (unsigned int i=0;i<8;++i)
335 {
336 kids[i]->doLeafWalk_const(c,v);
337 }
338 }
339 }
340
341 void OctCell::doLeafWalk(cellfn c, void* v)
342 {
343 if (leaf)
344 {
345 c(*this, v);
346 }
347 else
348 {
349 for (unsigned int i=0;i<8;++i)
350 {
351 kids[i]->doLeafWalk(c,v);
352 }
353 }
354 }
355
356 void OctCell::doLeafWalkWithKids_const(const_cellfn2 c, int k, void* v) const
357 {
358 if (leaf)
359 {
360 c(*this, k, v);
361 }
362 else
363 {
364 for (unsigned int i=0;i<8;++i)
365 {
366 kids[i]->doLeafWalkWithKids_const(c,i, v);
367 }
368 }
369
370
371
372 }
373
374
375 void OctCell::debug(bool fromroot)
376 {
377 linkCheck(fromroot);
378 }
379
380
381 void OctCell::linkCheck(bool fromroot)
382 {
383 if (fromroot && parent)
384 {
385 parent->debug(false);
386 }
387 else
388 {
389 if (leaf)
390 {
391 // if (leafinfo->owner!=this)
392 // {
393 // cout << this << " fails check\n";
394 // }
395 //cerr << centre[0] << ", " << centre[1] << ", " << centre[2] << "\n";
396 for (unsigned i=0;i<6;++i)
397 {
398 if (leafinfo->next[i])
399 {
400 //cerr << i << ',';
401 if (leafinfo->next[i]->owner->depth>10)
402 {
403 cout << this << " fails check ...\n";
404 }
405 }
406 }
407 //cerr << endl;
408
409 }
410 else
411 {
412 for (unsigned int i=0;i<8;++i)
413 {
414 kids[i]->debug(false);
415 }
416
417 }
418 }
419
420
421 }
422
423 bool OctCell::whohas(LeafInfo* li, bool fromroot)
424 {
425 if (fromroot && parent)
426 {
427 return parent->whohas(li, true);
428 }
429 else
430 {
431 bool has=false;
432 if (leaf)
433 {
434 if (leafinfo==li)
435 {
436 cerr << this << " has it as info\n";
437 has=true;
438 }
439 else
440 {
441 // if (leafinfo->owner!=this)
442 // {
443 // cout << this << " fails check\n";
444 // }
445 //cerr << centre[0] << ", " << centre[1] << ", " << centre[2] << "\n";
446 for (unsigned i=0;i<6;++i)
447 {
448 if (leafinfo->next[i]==li)
449 {
450 has=true;
451 cerr << this << " has it as " << i << " face. (d=";
452 cerr << depth << ") ";
453 cerr << centre[0] << ", " << centre[1] << ", " << centre[2] << endl;
454 }
455 }
456
457 }
458 //cerr << endl;
459
460 }
461 else
462 {
463 for (unsigned int i=0;i<8;++i)
464 {
465 has|=kids[i]->whohas(li, false);
466 }
467
468 }
469 return has;
470 }
471 }
472
473 namespace
474 {
475 void countleaves(const OctCell& c, void* v)
476 {
477 (*reinterpret_cast<int*>(v))++;
478 }
479
480
481 void writePoints(const OctCell& c, void* v)
482 {
483
484 ostream& os=*(reinterpret_cast<ostream*>(v));
485 static int i=1;
486 os << i++ << " " << (c.centre[0]-c.sides[0]/2) << ' ' << (c.centre[1]-c.sides[1]/2) << ' ' << (c.centre[2]-c.sides[2]/2) << '\n';
487 os << i++ << " " << (c.centre[0]+c.sides[0]/2) << ' ' << (c.centre[1]-c.sides[1]/2) << ' ' << (c.centre[2]-c.sides[2]/2) << '\n';
488 os << i++ << " " << (c.centre[0]+c.sides[0]/2) << ' ' << (c.centre[1]+c.sides[1]/2) << ' ' << (c.centre[2]-c.sides[2]/2) << '\n';
489 os << i++ << " " << (c.centre[0]-c.sides[0]/2) << ' ' << (c.centre[1]+c.sides[1]/2) << ' ' << (c.centre[2]-c.sides[2]/2) << '\n';
490
491 os << i++ << " " << (c.centre[0]-c.sides[0]/2) << ' ' << (c.centre[1]-c.sides[1]/2) << ' ' << (c.centre[2]+c.sides[2]/2) << '\n';
492 os << i++ << " " << (c.centre[0]+c.sides[0]/2) << ' ' << (c.centre[1]-c.sides[1]/2) << ' ' << (c.centre[2]+c.sides[2]/2) << '\n';
493 os << i++ << " " << (c.centre[0]+c.sides[0]/2) << ' ' << (c.centre[1]+c.sides[1]/2) << ' ' << (c.centre[2]+c.sides[2]/2) << '\n';
494 os << i++ << " " << (c.centre[0]-c.sides[0]/2) << ' ' << (c.centre[1]+c.sides[1]/2) << ' ' << (c.centre[2]+c.sides[2]/2) << '\n';
495 }
496
497 void writeEl(const OctCell& c, void* v)
498 {
499 static int e=1;
500 static int i=1;
501 ostream& os=reinterpret_cast<ostream&>(v);
502 os << e << " 5 0 "; // element number
503 os << i << ' ' << (i+1) << ' ' << (i+2) << ' ' << (i+3) << ' ';
504 os << (i+4) << ' ' << (i+5) << ' ' << (i+6) << ' ' << (i+7) << endl;
505 e++;
506 i+=8;
507 }
508
509
510 }
511
512 void OctCell::gmshDump()
513 {
514 if (parent)
515 {
516 parent->gmshDump();
517 }
518 else
519 {
520 ostream& os=cerr;
521 os << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n";
522
523 int c=0;
524 doLeafWalk_const(countleaves, &c);
525 os << c*8 << endl;
526 int i=1;
527 doLeafWalk_const(writePoints, &os);
528 os << "$EndNodes\n";
529 os << "$Elements\n";
530 os << c << endl;
531 doLeafWalk_const(writeEl, &os);
532 i=1;
533 os << "$EndElements\n";
534 }
535 }

  ViewVC Help
Powered by ViewVC 1.1.26