Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Stepniewski
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <cassert>
34 : #include "octahedron_rules.h"
35 : #include "rule_util.h"
36 : #include "grid_object_ids.h"
37 :
38 : namespace ug{
39 : namespace oct_rules
40 : {
41 :
42 : /// Output are the vertices of a octahedron rotated around its vertical axis
43 0 : void RotateOctahedron(int vrtsOut[NUM_VERTICES], int steps)
44 : {
45 0 : if(steps < 0)
46 0 : steps = (4 - ((-steps) % 4));
47 :
48 0 : for(int i = 0; i < 4; ++i)
49 0 : vrtsOut[((i + steps) % 4) + 1] = i+1;
50 :
51 0 : vrtsOut[0] = 0;
52 0 : vrtsOut[5] = 5;
53 0 : }
54 :
55 :
56 0 : int Refine(int* newIndsOut, int* newEdgeVrts, bool& newCenterOut, vector3* corners,
57 : bool*)
58 : {
59 0 : newCenterOut = false;
60 :
61 : // If a refinement rule is not implemented, fillCount will stay at 0.
62 : // Before returning, we'll check, whether fillCount is at 0 and will
63 : // perform refinement with an additional vertex in this case.
64 :
65 : // corner status is used to mark vertices, which are connected to refined edges
66 : // the value tells, how many associated edges are refined.
67 : int cornerStatus[6] = {0, 0, 0, 0, 0, 0};
68 :
69 : // here we'll store the index of each edge, which will be refined.
70 : // Size of the valid sub-array is numNewVrts, which is calculated below.
71 : // int refEdgeInds[NUM_EDGES];
72 :
73 : // count the number of new vertices and fill newVrtEdgeInds
74 : int numNewVrts = 0;
75 0 : for(int i = 0; i < NUM_EDGES; ++i){
76 0 : if(newEdgeVrts[i]){
77 : // refEdgeInds[numNewVrts] = i;
78 0 : ++numNewVrts;
79 :
80 : // adjust corner status of associated vertices
81 : const int* evi = EDGE_VRT_INDS[i];
82 : ++cornerStatus[evi[0]];
83 : ++cornerStatus[evi[1]];
84 : }
85 : }
86 :
87 : // the fillCount tells how much data has already been written to newIndsOut.
88 : int fillCount = 0;
89 :
90 :
91 : // convenience - indices where new edge-vrts, new face-vrts and new vol-vrts begin.
92 : const int E = NUM_VERTICES;
93 : const int V = NUM_VERTICES + NUM_EDGES + NUM_FACES;
94 :
95 : // depending on the number of new vertices, we will now apply different
96 : // refinement rules. Further distinction may have to be done.
97 0 : switch(numNewVrts){
98 0 : case 0:
99 : {
100 : // simply put the default octahedron back to newIndsOut
101 0 : newIndsOut[fillCount++] = GOID_OCTAHEDRON;
102 0 : newIndsOut[fillCount++] = 0;
103 0 : newIndsOut[fillCount++] = 1;
104 0 : newIndsOut[fillCount++] = 2;
105 0 : newIndsOut[fillCount++] = 3;
106 0 : newIndsOut[fillCount++] = 4;
107 0 : newIndsOut[fillCount++] = 5;
108 : }break;
109 :
110 : // TODO: Adaptive refinement cases still to be implemented!
111 :
112 : // REGULAR REFINEMENT
113 0 : case 12:
114 : {
115 : // We have to create 6 octahedrons and 8 tetrahedrons
116 : int& fi = fillCount;
117 : int* inds = newIndsOut;
118 :
119 : // lower tetrahedrons
120 0 : inds[fi++] = GOID_TETRAHEDRON;
121 0 : inds[fi++] = E + 4; inds[fi++] = E + 1;
122 0 : inds[fi++] = E; inds[fi++] = V;
123 :
124 0 : inds[fi++] = GOID_TETRAHEDRON;
125 0 : inds[fi++] = E + 5; inds[fi++] = E + 2;
126 0 : inds[fi++] = E + 1; inds[fi++] = V;
127 :
128 0 : inds[fi++] = GOID_TETRAHEDRON;
129 0 : inds[fi++] = E + 6; inds[fi++] = E + 3;
130 0 : inds[fi++] = E + 2; inds[fi++] = V;
131 :
132 0 : inds[fi++] = GOID_TETRAHEDRON;
133 0 : inds[fi++] = E + 7; inds[fi++] = E;
134 0 : inds[fi++] = E + 3; inds[fi++] = V;
135 :
136 : // upper tetrahedrons
137 0 : inds[fi++] = GOID_TETRAHEDRON;
138 0 : inds[fi++] = E + 9; inds[fi++] = E + 4;
139 0 : inds[fi++] = E + 8; inds[fi++] = V;
140 :
141 0 : inds[fi++] = GOID_TETRAHEDRON;
142 0 : inds[fi++] = E + 10; inds[fi++] = E + 5;
143 0 : inds[fi++] = E + 9; inds[fi++] = V;
144 :
145 0 : inds[fi++] = GOID_TETRAHEDRON;
146 0 : inds[fi++] = E + 11; inds[fi++] = E + 6;
147 0 : inds[fi++] = E + 10; inds[fi++] = V;
148 :
149 0 : inds[fi++] = GOID_TETRAHEDRON;
150 0 : inds[fi++] = E + 8; inds[fi++] = E + 7;
151 0 : inds[fi++] = E + 11; inds[fi++] = V;
152 :
153 : // For the 6 octahedrons we'll choose the shortest diagonals
154 : // and order the octahedral nodes, so that, the implicit
155 : // subdivision of the octahedron into tetrahedrons during
156 : // discretization happens alongside this shortest diagonal.
157 :
158 : int bestDiag[6];
159 :
160 : // for(int i = 0; i < 6; ++i)
161 : // bestDiag[i] = 2;
162 : //
163 : // if(corners){
164 : // //
165 : // // Calculate edge centers and element center
166 : // //
167 : //
168 : // vector3 c06, c07, c08, c09;
169 : // vector3 c10, c11, c12, c13;
170 : // vector3 c14, c15, c16, c17;
171 : // vector3 c18;
172 : //
173 : // // Bottom hemisphere edges
174 : // VecAdd(c06, corners[0], corners[1]);
175 : // VecScale(c06, c06, 0.5);
176 : //
177 : // VecAdd(c07, corners[0], corners[2]);
178 : // VecScale(c07, c07, 0.5);
179 : //
180 : // VecAdd(c08, corners[0], corners[3]);
181 : // VecScale(c08, c08, 0.5);
182 : //
183 : // VecAdd(c09, corners[0], corners[4]);
184 : // VecScale(c09, c09, 0.5);
185 : //
186 : // // Middle ring edges
187 : // VecAdd(c10, corners[1], corners[2]);
188 : // VecScale(c10, c10, 0.5);
189 : //
190 : // VecAdd(c11, corners[2], corners[3]);
191 : // VecScale(c11, c11, 0.5);
192 : //
193 : // VecAdd(c12, corners[3], corners[4]);
194 : // VecScale(c12, c12, 0.5);
195 : //
196 : // VecAdd(c13, corners[4], corners[1]);
197 : // VecScale(c13, c13, 0.5);
198 : //
199 : // // Top hemisphere edges
200 : // VecAdd(c14, corners[1], corners[5]);
201 : // VecScale(c14, c14, 0.5);
202 : //
203 : // VecAdd(c15, corners[2], corners[5]);
204 : // VecScale(c15, c15, 0.5);
205 : //
206 : // VecAdd(c16, corners[3], corners[5]);
207 : // VecScale(c16, c16, 0.5);
208 : //
209 : // VecAdd(c17, corners[4], corners[5]);
210 : // VecScale(c17, c17, 0.5);
211 : //
212 : // // Element center
213 : // VecAdd(c18, corners[0], corners[5]);
214 : // VecAppend(c18, corners[1], corners[2], corners[3], corners[4]);
215 : // VecScale(c18, c18, 1.0/6.0);
216 : //
217 : //
218 : // //
219 : // // Calculate the three diagonals of all 6 sub-octahedrons
220 : // //
221 : //
222 : // // Sub-octahedron 0
223 : // number oct0_d0_0608 = VecDistanceSq(c06, c08);
224 : // number oct0_d1_0709 = VecDistanceSq(c07, c09);
225 : // number oct0_d2_0018 = VecDistanceSq(corners[0], c18);
226 : //
227 : // number d = oct0_d2_0018;
228 : //
229 : // if(oct0_d1_0709 < d)
230 : // {
231 : // bestDiag[0] = 1;
232 : // d = oct0_d1_0709;
233 : // }
234 : //
235 : // if(oct0_d0_0608 < d)
236 : // {
237 : // bestDiag[0] = 0;
238 : // }
239 : //
240 : // //UG_LOG("oct0_d0_0608: " << oct0_d0_0608 << "; \t oct0_d1_0709: " << oct0_d1_0709 << "; \t oct0_d2_0018: " << oct0_d2_0018 << "; \t bestDiag0: " << bestDiag[0] << std::endl);
241 : //
242 : // // Sub-octahedron 1
243 : // number oct1_d0_0118 = VecDistanceSq(corners[1], c18);
244 : // number oct1_d1_1013 = VecDistanceSq(c10, c13);
245 : // number oct1_d2_0614 = VecDistanceSq(c06, c14);
246 : //
247 : // d = oct1_d2_0614;
248 : //
249 : // if(oct1_d1_1013 < d)
250 : // {
251 : // bestDiag[1] = 1;
252 : // d = oct1_d1_1013;
253 : // }
254 : //
255 : // if(oct1_d0_0118 < d)
256 : // {
257 : // bestDiag[1] = 0;
258 : // }
259 : //
260 : // //UG_LOG("oct1_d0_0118: " << oct1_d0_0118 << "; \t oct1_d1_1013: " << oct1_d1_1013 << "; \t oct1_d2_0614: " << oct1_d2_0614 << "; \t bestDiag1: " << bestDiag[1] << std::endl);
261 : //
262 : // // Sub-octahedron 2
263 : // number oct2_d0_1011 = VecDistanceSq(c10, c11);
264 : // number oct2_d1_0218 = VecDistanceSq(corners[2], c18);
265 : // number oct2_d2_0715 = VecDistanceSq(c07, c15);
266 : //
267 : // d = oct2_d2_0715;
268 : //
269 : // if(oct2_d1_0218 < d)
270 : // {
271 : // bestDiag[2] = 1;
272 : // d = oct2_d1_0218;
273 : // }
274 : //
275 : // if(oct2_d0_1011 < d)
276 : // {
277 : // bestDiag[2] = 0;
278 : // }
279 : //
280 : // //UG_LOG("oct2_d0_1011: " << oct2_d0_1011 << "; \t oct2_d1_0218: " << oct2_d1_0218 << "; \t oct2_d2_0715: " << oct2_d2_0715 << "; \t bestDiag2: " << bestDiag[2] << std::endl);
281 : //
282 : // // Sub-octahedron 3
283 : // number oct3_d0_1803 = VecDistanceSq(c18, corners[3]);
284 : // number oct3_d1_1112 = VecDistanceSq(c11, c12);
285 : // number oct3_d2_0816 = VecDistanceSq(c08, c16);
286 : //
287 : // d = oct3_d2_0816;
288 : //
289 : // if(oct3_d1_1112 < d)
290 : // {
291 : // bestDiag[3] = 1;
292 : // d = oct3_d1_1112;
293 : // }
294 : //
295 : // if(oct3_d0_1803 < d)
296 : // {
297 : // bestDiag[3] = 0;
298 : // }
299 : //
300 : // //UG_LOG("oct3_d0_1803: " << oct3_d0_1803 << "; \t oct3_d1_1112: " << oct3_d1_1112 << "; \t oct3_d2_0816: " << oct3_d2_0816 << "; \t bestDiag3: " << bestDiag[3] << std::endl);
301 : //
302 : // // Sub-octahedron 4
303 : // number oct4_d0_1312 = VecDistanceSq(c13, c12);
304 : // number oct4_d1_1804 = VecDistanceSq(c18, corners[4]);
305 : // number oct4_d2_0917 = VecDistanceSq(c09, c17);
306 : //
307 : // d = oct4_d2_0917;
308 : //
309 : // if(oct4_d1_1804 < d)
310 : // {
311 : // bestDiag[4] = 1;
312 : // d = oct4_d1_1804;
313 : // }
314 : //
315 : // if(oct4_d0_1312 < d)
316 : // {
317 : // bestDiag[4] = 0;
318 : // }
319 : //
320 : // //UG_LOG("oct4_d0_1312: " << oct4_d0_1312 << "; \t oct4_d1_1804: " << oct4_d1_1804 << "; \t oct4_d2_0917: " << oct4_d2_0917 << "; \t bestDiag4: " << bestDiag[4] << std::endl);
321 : //
322 : //
323 : // // Sub-octahedron 5
324 : // number oct5_d0_1416 = VecDistanceSq(c14, c16);
325 : // number oct5_d1_1517 = VecDistanceSq(c15, c17);
326 : // number oct5_d2_1805 = VecDistanceSq(c18, corners[5]);
327 : //
328 : // d = oct5_d0_1416;
329 : //
330 : // if(oct5_d1_1517 < d)
331 : // {
332 : // bestDiag[5] = 1;
333 : // d = oct5_d1_1517;
334 : // }
335 : //
336 : // if(oct5_d2_1805 < d)
337 : // {
338 : // bestDiag[5] = 0;
339 : // }
340 : //
341 : // //UG_LOG("oct5_d0_1416: " << oct5_d0_1416 << "; \t oct5_d1_1517: " << oct5_d1_1517 << "; \t oct5_d2_1805: " << oct5_d2_1805 << "; \t bestDiag5: " << bestDiag[5] << std::endl);
342 : // }
343 :
344 : // Fixed choice of the bestDiag (same as father's)
345 0 : for(int i = 0; i < 6; ++i)
346 0 : bestDiag[i] = 0;
347 :
348 : //
349 : // define octahedrons
350 : //
351 :
352 : // Sub-octahedron 0
353 0 : switch(bestDiag[0]){
354 0 : case 0: // diag 6-8
355 0 : inds[fi++] = GOID_OCTAHEDRON;
356 0 : inds[fi++] = 0; inds[fi++] = E;
357 0 : inds[fi++] = E + 1; inds[fi++] = E + 2;
358 0 : inds[fi++] = E + 3; inds[fi++] = V;
359 0 : break;
360 :
361 0 : case 1: // diag 7-9
362 0 : inds[fi++] = GOID_OCTAHEDRON;
363 0 : inds[fi++] = 0; inds[fi++] = E + 3;
364 0 : inds[fi++] = E; inds[fi++] = E + 1;
365 0 : inds[fi++] = E + 2; inds[fi++] = V;
366 0 : break;
367 :
368 0 : case 2: // diag 0-18
369 0 : inds[fi++] = GOID_OCTAHEDRON;
370 0 : inds[fi++] = E; inds[fi++] = V;
371 0 : inds[fi++] = E + 1; inds[fi++] = 0;
372 0 : inds[fi++] = E + 3; inds[fi++] = E + 2;
373 0 : break;
374 : }
375 :
376 : // Sub-octahedron 1
377 0 : switch(bestDiag[1]){
378 0 : case 0: // diag 1-18
379 0 : inds[fi++] = GOID_OCTAHEDRON;
380 0 : inds[fi++] = E; inds[fi++] = 1;
381 0 : inds[fi++] = E + 4; inds[fi++] = V;
382 0 : inds[fi++] = E + 7; inds[fi++] = E + 8;
383 0 : break;
384 :
385 0 : case 1: // diag 10-13
386 0 : inds[fi++] = GOID_OCTAHEDRON;
387 0 : inds[fi++] = E; inds[fi++] = E + 7;
388 0 : inds[fi++] = 1; inds[fi++] = E + 4;
389 0 : inds[fi++] = V; inds[fi++] = E + 8;
390 0 : break;
391 :
392 0 : case 2: // diag 6-14
393 0 : inds[fi++] = GOID_OCTAHEDRON;
394 0 : inds[fi++] = 1; inds[fi++] = E + 8;
395 0 : inds[fi++] = E + 4; inds[fi++] = E;
396 0 : inds[fi++] = E + 7; inds[fi++] = V;
397 0 : break;
398 : }
399 :
400 : // Sub-octahedron 2
401 0 : switch(bestDiag[2]){
402 0 : case 0: // diag 10-11
403 0 : inds[fi++] = GOID_OCTAHEDRON;
404 0 : inds[fi++] = E + 1; inds[fi++] = E + 4;
405 0 : inds[fi++] = 2; inds[fi++] = E + 5;
406 0 : inds[fi++] = V; inds[fi++] = E + 9;
407 0 : break;
408 :
409 0 : case 1: // diag 2-18
410 0 : inds[fi++] = GOID_OCTAHEDRON;
411 0 : inds[fi++] = E + 1; inds[fi++] = V;
412 0 : inds[fi++] = E + 4; inds[fi++] = 2;
413 0 : inds[fi++] = E + 5; inds[fi++] = E + 9;
414 0 : break;
415 :
416 0 : case 2: // diag 7-15
417 0 : inds[fi++] = GOID_OCTAHEDRON;
418 0 : inds[fi++] = E + 4; inds[fi++] = E + 9;
419 0 : inds[fi++] = 2; inds[fi++] = E + 1;
420 0 : inds[fi++] = V; inds[fi++] = E + 5;
421 0 : break;
422 : }
423 :
424 : // Sub-octahedron 3
425 0 : switch(bestDiag[3]){
426 0 : case 0: // diag 18-3
427 0 : inds[fi++] = GOID_OCTAHEDRON;
428 0 : inds[fi++] = E + 2; inds[fi++] = V;
429 0 : inds[fi++] = E + 5; inds[fi++] = 3;
430 0 : inds[fi++] = E + 6; inds[fi++] = E + 10;
431 0 : break;
432 :
433 0 : case 1: // diag 11-12
434 0 : inds[fi++] = GOID_OCTAHEDRON;
435 0 : inds[fi++] = E + 2; inds[fi++] = E + 6;
436 0 : inds[fi++] = V; inds[fi++] = E + 5;
437 0 : inds[fi++] = 3; inds[fi++] = E + 10;
438 0 : break;
439 :
440 0 : case 2: // diag 8-16
441 0 : inds[fi++] = GOID_OCTAHEDRON;
442 0 : inds[fi++] = V; inds[fi++] = E + 10;
443 0 : inds[fi++] = E + 5; inds[fi++] = E + 2;
444 0 : inds[fi++] = E + 6; inds[fi++] = 3;
445 0 : break;
446 : }
447 :
448 : // Sub-octahedron 4
449 0 : switch(bestDiag[4]){
450 0 : case 0: // diag 13-12
451 0 : inds[fi++] = GOID_OCTAHEDRON;
452 0 : inds[fi++] = E + 3; inds[fi++] = E + 7;
453 0 : inds[fi++] = V; inds[fi++] = E + 6;
454 0 : inds[fi++] = 4; inds[fi++] = E + 11;
455 0 : break;
456 :
457 0 : case 1: // diag 18-4
458 0 : inds[fi++] = GOID_OCTAHEDRON;
459 0 : inds[fi++] = E + 3; inds[fi++] = 4;
460 0 : inds[fi++] = E + 7; inds[fi++] = V;
461 0 : inds[fi++] = E + 6; inds[fi++] = E + 11;
462 0 : break;
463 :
464 0 : case 2: // diag 9-17
465 0 : inds[fi++] = GOID_OCTAHEDRON;
466 0 : inds[fi++] = E + 7; inds[fi++] = E + 11;
467 0 : inds[fi++] = V; inds[fi++] = E + 3;
468 0 : inds[fi++] = 4; inds[fi++] = E + 6;
469 0 : break;
470 : }
471 :
472 : // Sub-octahedron 5
473 0 : switch(bestDiag[5]){
474 0 : case 0: // diag 14-16
475 0 : inds[fi++] = GOID_OCTAHEDRON;
476 0 : inds[fi++] = V; inds[fi++] = E + 8;
477 0 : inds[fi++] = E + 9; inds[fi++] = E + 10;
478 0 : inds[fi++] = E + 11; inds[fi++] = 5;
479 0 : break;
480 :
481 0 : case 1: // diag 15-17
482 0 : inds[fi++] = GOID_OCTAHEDRON;
483 0 : inds[fi++] = V; inds[fi++] = E + 11;
484 0 : inds[fi++] = E + 8; inds[fi++] = E + 9;
485 0 : inds[fi++] = E + 10; inds[fi++] = 5;
486 0 : break;
487 :
488 0 : case 2: // diag 18-5
489 0 : inds[fi++] = GOID_OCTAHEDRON;
490 0 : inds[fi++] = E + 8; inds[fi++] = 5;
491 0 : inds[fi++] = E + 9; inds[fi++] = V;
492 0 : inds[fi++] = E + 11; inds[fi++] = E + 10;
493 0 : break;
494 : }
495 :
496 : // the rule requires a new center vertex
497 0 : newCenterOut = true;
498 :
499 : }break;
500 : }
501 :
502 : if(fillCount == 0){
503 : // call recursive refine
504 0 : fillCount = shared_rules::RecursiveRefine(newIndsOut, newEdgeVrts,
505 : FACE_VRT_INDS, FACE_EDGE_INDS, NUM_VERTICES,
506 : NUM_EDGES, NUM_FACES);
507 :
508 : // the rule requires a new center vertex
509 0 : newCenterOut = true;
510 : }
511 :
512 0 : return fillCount;
513 : }
514 :
515 :
516 0 : bool IsRegularRefRule(const int edgeMarks)
517 : {
518 : //todo: Think about regular octahedron ref rules.
519 0 : return false;
520 : }
521 :
522 : }// end of namespace oct_rules
523 : }// end of namespace ug
|
