Line data Source code
1 : /*
2 : * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
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 "tetrahedron_rules.h"
35 : #include "rule_util.h"
36 : #include "grid_object_ids.h"
37 :
38 : namespace ug{
39 : namespace tet_rules
40 : {
41 :
42 : /// global refinement rule information switching between regular and subdivision volume refinement
43 : static GlobalRefinementRule g_refinementRule = STANDARD;
44 :
45 0 : void SetRefinementRule(GlobalRefinementRule refRule)
46 : {
47 0 : g_refinementRule = refRule;
48 0 : }
49 :
50 0 : GlobalRefinementRule GetRefinementRule()
51 : {
52 0 : return g_refinementRule;
53 : }
54 :
55 :
56 : /// Rotates the given tetrahedron while keeping the specified point fixed
57 : /** The face opposed to the fixed point is rotated in counter-clockwise order when
58 : * viewed from the fixed point.
59 : * e.g.
60 : * \code
61 : * int vrts[] = {0, 1, 2, 3};
62 : * TetRotation(vrts, 3, 1);
63 : * // -> vrts == {2, 0, 1, 3}
64 : * \endcode
65 : */
66 0 : void TetRotation (
67 : int vrtsInOut[NUM_VERTICES],
68 : const int fixedPoint,
69 : const size_t steps)
70 : {
71 : // get the opposing face of the fixed point
72 0 : const int opFace = OPPOSED_OBJECT[fixedPoint][1];
73 0 : const int* finds = FACE_VRT_INDS[opFace];
74 :
75 0 : for(size_t i = 0; i < steps; ++i){
76 0 : const int tmp = vrtsInOut[finds[0]];
77 0 : vrtsInOut[finds[0]] = vrtsInOut[finds[2]];
78 0 : vrtsInOut[finds[2]] = vrtsInOut[finds[1]];
79 0 : vrtsInOut[finds[1]] = tmp;
80 : }
81 0 : }
82 :
83 0 : void InverseTetTransform(int* indsOut, const int* transformedInds){
84 : UG_ASSERT(indsOut != transformedInds, "The arrays have to differ!");
85 0 : for(int i = 0; i < NUM_VERTICES; ++i)
86 0 : indsOut[transformedInds[i]] = i;
87 0 : }
88 :
89 :
90 0 : int Refine(int* newIndsOut, int* newEdgeVrts, bool& newCenterOut,
91 : vector3* corners, bool*)
92 : {
93 0 : newCenterOut = false;
94 : // If a refinement rule is not implemented, fillCount will stay at 0.
95 : // Before returning, we'll check, whether fillCount is at 0 and will
96 : // perform refinement with an additional vertex in this case.
97 :
98 : // corner status is used to mark vertices, which are connected to refined edges
99 : // the value tells, how many associated edges are refined.
100 0 : int cornerStatus[4] = {0, 0, 0, 0};
101 :
102 : // here we'll store the index of each edge, which will be refined.
103 : // Size of the valid sub-array is numNewVrts, which is calculated below.
104 : int refEdgeInds[NUM_EDGES];
105 :
106 : // count the number of new vertices and fill newVrtEdgeInds
107 : int numNewVrts = 0;
108 0 : for(int i = 0; i < NUM_EDGES; ++i){
109 0 : if(newEdgeVrts[i]){
110 0 : refEdgeInds[numNewVrts] = i;
111 0 : ++numNewVrts;
112 :
113 : // adjust corner status of associated vertices
114 0 : const int* evi = EDGE_VRT_INDS[i];
115 0 : ++cornerStatus[evi[0]];
116 0 : ++cornerStatus[evi[1]];
117 : }
118 : }
119 :
120 : // the fillCount tells how much data has already been written to newIndsOut.
121 : int fillCount = 0;
122 :
123 : // depending on the number of new vertices, we will now apply different
124 : // refinement rules. Further distinction may have to be done.
125 0 : switch(numNewVrts){
126 0 : case 0:
127 : {
128 : // simply put the default tetrahedron back to newIndsOut
129 0 : newIndsOut[fillCount++] = GOID_TETRAHEDRON;
130 0 : newIndsOut[fillCount++] = 0;
131 0 : newIndsOut[fillCount++] = 1;
132 0 : newIndsOut[fillCount++] = 2;
133 0 : newIndsOut[fillCount++] = 3;
134 : }break;
135 :
136 0 : case 1:
137 : {
138 0 : int refEdgeInd = refEdgeInds[0];
139 : // the two faces which are not connected to the refined edge
140 : // serve as bottom for the new tetrahedrons.
141 0 : for(int i_face = 0; i_face < NUM_FACES; ++i_face){
142 0 : if(!FACE_CONTAINS_EDGE[i_face][refEdgeInd]){
143 0 : const int* fvi = FACE_VRT_INDS[i_face];
144 :
145 0 : newIndsOut[fillCount++] = GOID_TETRAHEDRON;
146 0 : newIndsOut[fillCount++] = fvi[0];
147 0 : newIndsOut[fillCount++] = fvi[1];
148 0 : newIndsOut[fillCount++] = fvi[2];
149 0 : newIndsOut[fillCount++] = NUM_VERTICES + refEdgeInd;
150 : }
151 : }
152 : }break;
153 :
154 0 : case 2:
155 : {
156 : // two types exist: The two refined edges share a vertex or not.
157 0 : if(OPPOSED_EDGE[refEdgeInds[0]] == refEdgeInds[1]){
158 : // they do not share an edge
159 : // we create a local order from refEdgeInds[0], which
160 : // always has to be an edge in the base-triangle and
161 : // refEdgeInds[1], which always connects a base-corner with
162 : // the top.
163 0 : const int v0 = EDGE_VRT_INDS[refEdgeInds[0]][0];
164 0 : const int v1 = EDGE_VRT_INDS[refEdgeInds[0]][1];
165 0 : const int v2 = EDGE_VRT_INDS[refEdgeInds[1]][0];
166 0 : const int v3 = EDGE_VRT_INDS[refEdgeInds[1]][1];
167 0 : const int n0 = refEdgeInds[0] + NUM_VERTICES;
168 0 : const int n1 = refEdgeInds[1] + NUM_VERTICES;
169 :
170 : // from this local order we can now construct the 4 new tetrahedrons.
171 : int& fi = fillCount;
172 : int* inds = newIndsOut;
173 0 : inds[fi++] = GOID_TETRAHEDRON;
174 0 : inds[fi++] = v0; inds[fi++] = n0;
175 0 : inds[fi++] = v2; inds[fi++] = n1;
176 :
177 0 : inds[fi++] = GOID_TETRAHEDRON;
178 0 : inds[fi++] = v0; inds[fi++] = n0;
179 0 : inds[fi++] = n1; inds[fi++] = v3;
180 :
181 0 : inds[fi++] = GOID_TETRAHEDRON;
182 0 : inds[fi++] = n0; inds[fi++] = v1;
183 0 : inds[fi++] = v2; inds[fi++] = n1;
184 :
185 0 : inds[fi++] = GOID_TETRAHEDRON;
186 0 : inds[fi++] = n0; inds[fi++] = v1;
187 0 : inds[fi++] = n1; inds[fi++] = v3;
188 : }
189 : else{
190 : // they share an edge
191 : // We have to create a pyramid and a tetrahedron.
192 : // get the triangle, which contains both refEdges
193 0 : int refTriInd = FACE_FROM_EDGES[refEdgeInds[0]][refEdgeInds[1]];
194 0 : const int* f = FACE_VRT_INDS[refTriInd];
195 :
196 : // find the edge (v0, v1) in refTri, which won't be refined
197 : int v0 = -1; int v1 = -1; int v2 = -1;
198 0 : for(int i = 0; i < 3; ++i){
199 0 : v0 = f[i];
200 0 : v1 = f[(i+1)%3];
201 0 : v2 = f[(i+2)%3];
202 0 : if(cornerStatus[v0] == 1 && cornerStatus[v1] == 1)
203 : break;
204 : }
205 :
206 : // make sure that v2 is connected to two refined edges
207 : assert(cornerStatus[v2] == 2);
208 :
209 : // get the new vertex on edge v1v2 and on v2v0
210 0 : int v1v2 = EDGE_FROM_VRTS[v1][v2] + NUM_VERTICES;
211 0 : int v2v0 = EDGE_FROM_VRTS[v2][v0] + NUM_VERTICES;
212 :
213 : // get the top (vertex with cornerState 0)
214 : int vtop = -1;
215 0 : for(int i = 0; i < NUM_VERTICES; ++i){
216 0 : if(cornerStatus[i] == 0){
217 : vtop = i;
218 : break;
219 : }
220 : }
221 :
222 : // now lets build the pyramid
223 : int& fi = fillCount;
224 : int* inds = newIndsOut;
225 0 : inds[fi++] = GOID_PYRAMID;
226 0 : inds[fi++] = v0; inds[fi++] = v1;
227 0 : inds[fi++] = v1v2; inds[fi++] = v2v0;
228 0 : inds[fi++] = vtop;
229 :
230 : // and now the terahedron
231 0 : inds[fi++] = GOID_TETRAHEDRON;
232 0 : inds[fi++] = v2; inds[fi++] = vtop;
233 0 : inds[fi++] = v2v0; inds[fi++] = v1v2;
234 : }
235 : }break;
236 :
237 0 : case 3:
238 : {
239 : // different possibilities exist. First we'll treat the one,
240 : // where all new vertices lie on the edges of one triangle.
241 : // Note that refTriInd could be -1 after the call.
242 0 : int refTriInd = FACE_FROM_EDGES[refEdgeInds[0]][refEdgeInds[1]];
243 0 : if(refTriInd == FACE_FROM_EDGES[refEdgeInds[1]][refEdgeInds[2]])
244 : {
245 : // all three lie in one plane (refTriInd has to be != -1 to get here)
246 : // get the top (vertex with cornerState 0)
247 : int vtop = -1;
248 0 : for(int i = 0; i < NUM_VERTICES; ++i){
249 0 : if(cornerStatus[i] == 0){
250 : vtop = i;
251 : break;
252 : }
253 : }
254 :
255 0 : const int* f = FACE_VRT_INDS[refTriInd];
256 :
257 : // create four new tetrahedrons
258 0 : const int v0 = f[0]; const int v1 = f[1]; const int v2 = f[2];
259 0 : const int v0v1 = EDGE_FROM_VRTS[v0][v1] + NUM_VERTICES;
260 0 : const int v1v2 = EDGE_FROM_VRTS[v1][v2] + NUM_VERTICES;
261 0 : const int v2v0 = EDGE_FROM_VRTS[v2][v0] + NUM_VERTICES;
262 :
263 : int& fi = fillCount;
264 : int* inds = newIndsOut;
265 0 : inds[fi++] = GOID_TETRAHEDRON;
266 0 : inds[fi++] = v0; inds[fi++] = vtop;
267 0 : inds[fi++] = v0v1; inds[fi++] = v2v0;
268 0 : inds[fi++] = GOID_TETRAHEDRON;
269 0 : inds[fi++] = v1; inds[fi++] = vtop;
270 0 : inds[fi++] = v1v2; inds[fi++] = v0v1;
271 0 : inds[fi++] = GOID_TETRAHEDRON;
272 0 : inds[fi++] = v2; inds[fi++] = vtop;
273 0 : inds[fi++] = v2v0; inds[fi++] = v1v2;
274 0 : inds[fi++] = GOID_TETRAHEDRON;
275 0 : inds[fi++] = v0v1; inds[fi++] = vtop;
276 0 : inds[fi++] = v1v2; inds[fi++] = v2v0;
277 : }
278 : else{
279 : // we have to further distinguish.
280 : // gather the corners with corner-state 3 and corner state 1
281 : int corner3 = -1;
282 : int freeCorner[NUM_VERTICES];
283 : int freeCount = 0;
284 0 : for(int i = 0; i < NUM_VERTICES; ++i){
285 0 : if(cornerStatus[i] == 3)
286 : corner3 = i;
287 : else
288 0 : freeCorner[freeCount++] = i;
289 : }
290 :
291 0 : if(corner3 != -1){
292 : // a corner with corner state 3 exists.
293 : assert(freeCount == 3);
294 :
295 : // get the face which won't be refined (required for correct orientation)
296 0 : int freeTriInd = FACE_FROM_VRTS[freeCorner[0]][freeCorner[1]]
297 0 : [freeCorner[2]];
298 :
299 : // the free tri is the new base, corner3 the top
300 0 : const int* f = FACE_VRT_INDS[freeTriInd];
301 0 : int v2v3 = EDGE_FROM_VRTS[f[2]][corner3] + NUM_VERTICES;
302 0 : int v1v3 = EDGE_FROM_VRTS[f[1]][corner3] + NUM_VERTICES;
303 0 : int v0v3 = EDGE_FROM_VRTS[f[0]][corner3] + NUM_VERTICES;
304 :
305 : // add a prism and a tetrahedron
306 : int& fi = fillCount;
307 : int* inds = newIndsOut;
308 0 : inds[fi++] = GOID_PRISM;
309 0 : inds[fi++] = f[0]; inds[fi++] = f[1]; inds[fi++] = f[2];
310 0 : inds[fi++] = v0v3; inds[fi++] = v1v3; inds[fi++] = v2v3;
311 :
312 0 : inds[fi++] = GOID_TETRAHEDRON;
313 0 : inds[fi++] = v2v3; inds[fi++] = corner3; inds[fi++] = v0v3;
314 0 : inds[fi++] = v1v3;
315 : }
316 : }
317 : }break;
318 :
319 : case 4:
320 : {
321 : // multiple settings with 4 refined edges exist.
322 : // Either two opposing edges are not marked (case all2 == true)
323 : // or the two unmarked edges are contained in 1 triangle
324 :
325 : // check whether all vertices have corner state 2
326 : bool all2 = true;
327 0 : for(int i = 0; i < NUM_VERTICES; ++i){
328 0 : if(cornerStatus[i] !=2){
329 : all2 = false;
330 : break;
331 : }
332 : }
333 :
334 0 : if(all2){
335 : // we've got a straight cut.
336 : // we'll rotate the tetrahedron around the tip so, that edge 2 won't be refined.
337 : int steps = 0;
338 0 : if(!newEdgeVrts[EDGE_FROM_VRTS[0][1]])
339 : steps = 2;
340 0 : else if(!newEdgeVrts[EDGE_FROM_VRTS[1][2]])
341 : steps = 1;
342 :
343 0 : int t[NUM_VERTICES] = {0, 1, 2, 3};
344 0 : TetRotation(t, 3, steps);
345 :
346 0 : const int v0v1 = EDGE_FROM_VRTS[t[0]][t[1]] + NUM_VERTICES;
347 0 : const int v1v2 = EDGE_FROM_VRTS[t[1]][t[2]] + NUM_VERTICES;
348 0 : const int v0v3 = EDGE_FROM_VRTS[t[0]][t[3]] + NUM_VERTICES;
349 0 : const int v2v3 = EDGE_FROM_VRTS[t[2]][t[3]] + NUM_VERTICES;
350 :
351 : assert(newEdgeVrts[v0v1 - NUM_VERTICES]);
352 : assert(newEdgeVrts[v1v2 - NUM_VERTICES]);
353 : assert(newEdgeVrts[v0v3 - NUM_VERTICES]);
354 : assert(newEdgeVrts[v2v3 - NUM_VERTICES]);
355 :
356 : // now build two prisms
357 : int& fi = fillCount;
358 : int* inds = newIndsOut;
359 :
360 0 : inds[fi++] = GOID_PRISM;
361 0 : inds[fi++] = t[0]; inds[fi++] = v0v3; inds[fi++] = v0v1;
362 0 : inds[fi++] = t[2]; inds[fi++] = v2v3; inds[fi++] = v1v2;
363 :
364 0 : inds[fi++] = GOID_PRISM;
365 0 : inds[fi++] = t[1]; inds[fi++] = v1v2; inds[fi++] = v0v1;
366 0 : inds[fi++] = t[3]; inds[fi++] = v2v3; inds[fi++] = v0v3;
367 : }
368 : else{
369 : // one corner has state 1, one has state 3 and two have state 2.
370 : // Rotate the tet so that 1 is at the top and 4 is at the origin
371 0 : int I[NUM_VERTICES] = {0, 1, 2, 3};
372 :
373 : int s1 = -1;
374 0 : for(int i = 0; i < 4; ++i){
375 0 : if(cornerStatus[i] == 1){
376 : s1 = i;
377 : break;
378 : }
379 : }
380 :
381 0 : if(s1 != 3){
382 0 : const int fixedPoint = (s1 + 2) % 3;
383 0 : TetRotation(I, fixedPoint, 1);
384 : }
385 :
386 0 : if(cornerStatus[I[1]] == 3)
387 0 : TetRotation(I, 3, 2);
388 0 : else if(cornerStatus[I[2]] == 3)
389 0 : TetRotation(I, 3, 1);
390 :
391 : // indices of edge-center vertices
392 0 : const int v01 = EDGE_FROM_VRTS[I[0]][I[1]] + NUM_VERTICES;
393 0 : const int v02 = EDGE_FROM_VRTS[I[0]][I[2]] + NUM_VERTICES;
394 0 : const int v03 = EDGE_FROM_VRTS[I[0]][I[3]] + NUM_VERTICES;
395 0 : const int v12 = EDGE_FROM_VRTS[I[1]][I[2]] + NUM_VERTICES;
396 :
397 : int& fi = fillCount;
398 : int* inds = newIndsOut;
399 :
400 0 : inds[fi++] = GOID_TETRAHEDRON;
401 0 : inds[fi++] = I[0]; inds[fi++] = v01;
402 0 : inds[fi++] = v02; inds[fi++] = v03;
403 :
404 0 : inds[fi++] = GOID_TETRAHEDRON;
405 0 : inds[fi++] = v01; inds[fi++] = v12;
406 0 : inds[fi++] = v02; inds[fi++] = v03;
407 :
408 0 : inds[fi++] = GOID_PYRAMID;
409 0 : inds[fi++] = I[3]; inds[fi++] = I[1];
410 0 : inds[fi++] = v01; inds[fi++] = v03;
411 0 : inds[fi++] = v12;
412 :
413 0 : inds[fi++] = GOID_PYRAMID;
414 0 : inds[fi++] = I[2]; inds[fi++] = I[3];
415 0 : inds[fi++] = v03; inds[fi++] = v02;
416 0 : inds[fi++] = v12;
417 : }
418 : }break;
419 :
420 : case 5:{
421 : // only one edge is not marked for refinement
422 : int unmarkedEdge = 0;
423 0 : for(int i = 0; i < NUM_EDGES; ++i){
424 0 : if(!newEdgeVrts[i]){
425 : unmarkedEdge = i;
426 : break;
427 : }
428 : }
429 :
430 0 : int uei0 = EDGE_VRT_INDS[unmarkedEdge][0];
431 0 : int uei1 = EDGE_VRT_INDS[unmarkedEdge][1];
432 :
433 : // orientate the tetrahedron so that the unmarked edge connects
434 : // vertices 2 and 3
435 0 : int I[NUM_VERTICES] = {0, 1, 2, 3};
436 : // if the unmarked edge lies in the base triangle, we'll first have to
437 : // rotate so that it connects a base-vertex with the top
438 0 : if(unmarkedEdge < 3){
439 0 : const int fixedPoint = (uei1 + 1) % 3;
440 0 : TetRotation(I, fixedPoint, 1);
441 : int IInv[4];
442 0 : InverseTetTransform(IInv, I);
443 :
444 0 : uei0 = IInv[uei0];
445 0 : uei1 = IInv[uei1];
446 0 : unmarkedEdge = EDGE_FROM_VRTS[uei0][uei1];
447 : }
448 :
449 : // now rotate the tet to its final place
450 0 : if(unmarkedEdge == 3)
451 0 : TetRotation(I, 3, 2);
452 0 : else if(unmarkedEdge == 4)
453 0 : TetRotation(I, 3, 1);
454 :
455 : // We obtained the final permutation I. Now create new elements
456 : // indices of edge-center vertices
457 0 : const int v01 = EDGE_FROM_VRTS[I[0]][I[1]] + NUM_VERTICES;
458 0 : const int v02 = EDGE_FROM_VRTS[I[0]][I[2]] + NUM_VERTICES;
459 0 : const int v03 = EDGE_FROM_VRTS[I[0]][I[3]] + NUM_VERTICES;
460 0 : const int v12 = EDGE_FROM_VRTS[I[1]][I[2]] + NUM_VERTICES;
461 0 : const int v13 = EDGE_FROM_VRTS[I[1]][I[3]] + NUM_VERTICES;
462 :
463 : int& fi = fillCount;
464 : int* inds = newIndsOut;
465 :
466 0 : inds[fi++] = GOID_TETRAHEDRON;
467 0 : inds[fi++] = I[0]; inds[fi++] = v01;
468 0 : inds[fi++] = v02; inds[fi++] = v03;
469 :
470 0 : inds[fi++] = GOID_TETRAHEDRON;
471 0 : inds[fi++] = I[1]; inds[fi++] = v12;
472 0 : inds[fi++] = v01; inds[fi++] = v13;
473 :
474 0 : inds[fi++] = GOID_PYRAMID;
475 0 : inds[fi++] = v02; inds[fi++] = v12;
476 0 : inds[fi++] = v13; inds[fi++] = v03;
477 0 : inds[fi++] = v01;
478 :
479 0 : inds[fi++] = GOID_PRISM;
480 0 : inds[fi++] = v02; inds[fi++] = v12; inds[fi++] = I[2];
481 0 : inds[fi++] = v03; inds[fi++] = v13; inds[fi++] = I[3];
482 :
483 :
484 : }break;
485 :
486 : // REGULAR REFINEMENT
487 0 : case 6:
488 : {
489 : // depending on the user specified global refinement rule, we will now apply different
490 : // refinement rules.
491 0 : switch(g_refinementRule){
492 0 : case STANDARD:
493 : {
494 : int& fi = fillCount;
495 : int* inds = newIndsOut;
496 :
497 0 : inds[fi++] = GOID_TETRAHEDRON;
498 0 : inds[fi++] = 0; inds[fi++] = NUM_VERTICES + 0;
499 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 3;
500 :
501 0 : inds[fi++] = GOID_TETRAHEDRON;
502 0 : inds[fi++] = 1; inds[fi++] = NUM_VERTICES + 1;
503 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 4;
504 :
505 0 : inds[fi++] = GOID_TETRAHEDRON;
506 0 : inds[fi++] = 2; inds[fi++] = NUM_VERTICES + 2;
507 0 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 5;
508 :
509 0 : inds[fi++] = GOID_TETRAHEDRON;
510 0 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 4;
511 0 : inds[fi++] = NUM_VERTICES + 5; inds[fi++] = 3;
512 :
513 : // for the remaining four tetrahedrons, we'll choose the shortest diagonal
514 : int bestDiag = 2;
515 0 : if(corners){
516 : // there are three diagonals between the following edge-centers:
517 : // 0-5, 1-3, 2-4
518 : vector3 c1, c2;
519 :
520 : // 0-5
521 : VecAdd(c1, corners[0], corners[1]);
522 : VecScale(c1, c1, 0.5);
523 : VecAdd(c2, corners[2], corners[3]);
524 : VecScale(c2, c2, 0.5);
525 0 : number d05 = VecDistanceSq(c1, c2);
526 :
527 : // 1-3
528 : VecAdd(c1, corners[1], corners[2]);
529 : VecScale(c1, c1, 0.5);
530 : VecAdd(c2, corners[0], corners[3]);
531 : VecScale(c2, c2, 0.5);
532 0 : number d13 = VecDistanceSq(c1, c2);
533 :
534 : // 2-4
535 : VecAdd(c1, corners[0], corners[2]);
536 : VecScale(c1, c1, 0.5);
537 : VecAdd(c2, corners[1], corners[3]);
538 : VecScale(c2, c2, 0.5);
539 0 : number d = VecDistanceSq(c1, c2);
540 :
541 0 : if(d13 < d){
542 : bestDiag = 1;
543 : d = d13;
544 : }
545 0 : if(d05 < d){
546 : bestDiag = 0;
547 : }
548 : }
549 :
550 0 : switch(bestDiag){
551 : case 0:// diag: 0-5
552 0 : inds[fi++] = GOID_TETRAHEDRON;
553 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 1;
554 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
555 :
556 0 : inds[fi++] = GOID_TETRAHEDRON;
557 0 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 4;
558 0 : inds[fi++] = NUM_VERTICES + 5; inds[fi++] = NUM_VERTICES + 0;
559 :
560 0 : inds[fi++] = GOID_TETRAHEDRON;
561 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
562 0 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 0;
563 :
564 0 : inds[fi++] = GOID_TETRAHEDRON;
565 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 3;
566 0 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 5;
567 :
568 : break;
569 :
570 0 : case 1:// diag: 1-3
571 0 : inds[fi++] = GOID_TETRAHEDRON;
572 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 1;
573 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 3;
574 :
575 0 : inds[fi++] = GOID_TETRAHEDRON;
576 0 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 4;
577 0 : inds[fi++] = NUM_VERTICES + 5; inds[fi++] = NUM_VERTICES + 3;
578 :
579 0 : inds[fi++] = GOID_TETRAHEDRON;
580 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
581 0 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 1;
582 :
583 0 : inds[fi++] = GOID_TETRAHEDRON;
584 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 3;
585 0 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 1;
586 :
587 : break;
588 :
589 0 : case 2:// diag 2-4
590 0 : inds[fi++] = GOID_TETRAHEDRON;
591 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 2;
592 0 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 4;
593 :
594 0 : inds[fi++] = GOID_TETRAHEDRON;
595 0 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 2;
596 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 4;
597 :
598 0 : inds[fi++] = GOID_TETRAHEDRON;
599 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 3;
600 0 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 5;
601 :
602 0 : inds[fi++] = GOID_TETRAHEDRON;
603 0 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 1;
604 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
605 :
606 : break;
607 : }
608 : }break;
609 :
610 0 : case HYBRID_TET_OCT:
611 : {
612 : // REGULAR REFINEMENT
613 : int& fi = fillCount;
614 : int* inds = newIndsOut;
615 :
616 : // outer tetrahedrons analogously defined to STANDARD case
617 0 : inds[fi++] = GOID_TETRAHEDRON;
618 0 : inds[fi++] = 0; inds[fi++] = NUM_VERTICES + 0;
619 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 3;
620 :
621 0 : inds[fi++] = GOID_TETRAHEDRON;
622 0 : inds[fi++] = 1; inds[fi++] = NUM_VERTICES + 1;
623 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 4;
624 :
625 0 : inds[fi++] = GOID_TETRAHEDRON;
626 0 : inds[fi++] = 2; inds[fi++] = NUM_VERTICES + 2;
627 0 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 5;
628 :
629 0 : inds[fi++] = GOID_TETRAHEDRON;
630 0 : inds[fi++] = NUM_VERTICES + 3; inds[fi++] = NUM_VERTICES + 4;
631 0 : inds[fi++] = NUM_VERTICES + 5; inds[fi++] = 3;
632 :
633 : // inner octahedron:
634 : // For the remaining inner cavity we'll choose the shortest diagonal
635 : // and order the octahedral nodes, so that, the implicit
636 : // subdivision of the octahedron into tetrahedrons during
637 : // discretization happens alongside exactly this shortest diagonal.
638 : int bestDiag = 2;
639 0 : if(corners){
640 : // there are three diagonals between the following edge-centers:
641 : // 0-5, 1-3, 2-4
642 : vector3 c1, c2;
643 :
644 : // 0-5
645 : VecAdd(c1, corners[0], corners[1]);
646 : VecScale(c1, c1, 0.5);
647 : VecAdd(c2, corners[2], corners[3]);
648 : VecScale(c2, c2, 0.5);
649 0 : number d05 = VecDistanceSq(c1, c2);
650 :
651 : // 1-3
652 : VecAdd(c1, corners[1], corners[2]);
653 : VecScale(c1, c1, 0.5);
654 : VecAdd(c2, corners[0], corners[3]);
655 : VecScale(c2, c2, 0.5);
656 0 : number d13 = VecDistanceSq(c1, c2);
657 :
658 : // 2-4
659 : VecAdd(c1, corners[0], corners[2]);
660 : VecScale(c1, c1, 0.5);
661 : VecAdd(c2, corners[1], corners[3]);
662 : VecScale(c2, c2, 0.5);
663 0 : number d = VecDistanceSq(c1, c2);
664 :
665 0 : if(d13 < d){
666 : bestDiag = 1;
667 : d = d13;
668 : }
669 0 : if(d05 < d){
670 : bestDiag = 0;
671 : }
672 : }
673 :
674 0 : switch(bestDiag){
675 : case 0:// diag: 0-5
676 0 : inds[fi++] = GOID_OCTAHEDRON;
677 0 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 0;
678 0 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 5;
679 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 3;
680 : break;
681 :
682 0 : case 1:// diag: 1-3
683 0 : inds[fi++] = GOID_OCTAHEDRON;
684 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 3;
685 0 : inds[fi++] = NUM_VERTICES + 4; inds[fi++] = NUM_VERTICES + 1;
686 0 : inds[fi++] = NUM_VERTICES + 2; inds[fi++] = NUM_VERTICES + 5;
687 : break;
688 :
689 0 : case 2:// diag 2-4
690 0 : inds[fi++] = GOID_OCTAHEDRON;
691 0 : inds[fi++] = NUM_VERTICES + 1; inds[fi++] = NUM_VERTICES + 2;
692 0 : inds[fi++] = NUM_VERTICES + 0; inds[fi++] = NUM_VERTICES + 4;
693 0 : inds[fi++] = NUM_VERTICES + 5; inds[fi++] = NUM_VERTICES + 3;
694 : break;
695 : }break;
696 :
697 : }break;
698 : }
699 : }break;
700 : }
701 :
702 0 : if(fillCount == 0){
703 : // call recursive refine
704 0 : fillCount = shared_rules::RecursiveRefine(newIndsOut, newEdgeVrts,
705 : FACE_VRT_INDS, FACE_EDGE_INDS, NUM_VERTICES,
706 : NUM_EDGES, NUM_FACES);
707 :
708 : // the rule requires a new center vertex
709 0 : newCenterOut = true;
710 : }
711 :
712 0 : return fillCount;
713 : }
714 :
715 :
716 0 : bool IsRegularRefRule(const int edgeMarks)
717 : {
718 0 : return false;
719 : }
720 :
721 : }// end of namespace tet_rules
722 : }// end of namespace ug
|
