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 "pyramid_rules.h"
35 : #include "rule_util.h"
36 : #include "grid_object_ids.h"
37 :
38 : namespace ug{
39 : namespace pyra_rules{
40 :
41 : /// Output are the vertices of a pyramid rotated around its vertical axis
42 0 : void RotatePyramid(int vrtsOut[NUM_VERTICES], int steps)
43 : {
44 0 : if(steps < 0)
45 0 : steps = (4 - ((-steps) % 4));
46 :
47 0 : for(int i = 0; i < 4; ++i)
48 0 : vrtsOut[(i + steps) % 4] = i;
49 :
50 0 : vrtsOut[4] = 4;
51 0 : }
52 :
53 :
54 :
55 :
56 0 : int Refine(int* newIndsOut, int* newEdgeVrts, bool& newCenterOut, vector3*, bool* isSnapPoint)
57 : {
58 0 : newCenterOut = false;
59 : // If a refinement rule is not implemented, fillCount will stay at 0.
60 : // Before returning, we'll check, whether fillCount is at 0 and will
61 : // perform refinement with an additional vertex in this case.
62 :
63 : // corner status is used to mark vertices, which are connected to refined edges
64 : // the value tells, how many associated edges are refined.
65 0 : int cornerStatus[5] = {0, 0, 0, 0, 0};
66 :
67 : // here we'll store the index of each edge, which will be refined.
68 : // Size of the valid sub-array is numNewVrts, which is calculated below.
69 : int refEdgeInds[NUM_EDGES];
70 :
71 : // count the number of new vertices and fill newVrtEdgeInds
72 : int numNewVrts = 0;
73 0 : for(int i = 0; i < NUM_EDGES; ++i){
74 0 : if(newEdgeVrts[i]){
75 0 : refEdgeInds[numNewVrts] = i;
76 0 : ++numNewVrts;
77 :
78 : // adjust corner status of associated vertices
79 0 : const int* evi = EDGE_VRT_INDS[i];
80 0 : ++cornerStatus[evi[0]];
81 0 : ++cornerStatus[evi[1]];
82 : }
83 : }
84 :
85 : // snap-point handling
86 : int numSnapPoints = 0;
87 0 : if(isSnapPoint){
88 0 : for(int i = 0; i < NUM_VERTICES; ++i){
89 0 : if(isSnapPoint[i])
90 0 : ++numSnapPoints;
91 : }
92 : }
93 :
94 : bool snapPointsProcessed = numSnapPoints == 0 ? true : false;
95 :
96 : // the fillCount tells how much data has already been written to newIndsOut.
97 : int fillCount = 0;
98 :
99 : // convenience - indices where new edge-vrts, new face-vrts and new vol-vrts begin.
100 : const int E = NUM_VERTICES;
101 : const int F = NUM_VERTICES + NUM_EDGES;
102 : // const int V = NUM_VERTICES + NUM_EDGES + NUM_FACES;
103 :
104 : // depending on the number of new vertices, we will now apply different
105 : // refinement rules. Further distinction may have to be done.
106 0 : switch(numNewVrts){
107 0 : case 0:
108 : {
109 : // simply put the default pyramid back to newIndsOut
110 0 : newIndsOut[fillCount++] = GOID_PYRAMID;
111 0 : newIndsOut[fillCount++] = 0;
112 0 : newIndsOut[fillCount++] = 1;
113 0 : newIndsOut[fillCount++] = 2;
114 0 : newIndsOut[fillCount++] = 3;
115 0 : newIndsOut[fillCount++] = 4;
116 0 : }break;
117 :
118 0 : case 1: // only one edge has to be refined.
119 : {
120 0 : const int refEdge = refEdgeInds[0];
121 0 : const int nVrt = refEdge + E;
122 :
123 : int& fi = fillCount;
124 : int* inds = newIndsOut;
125 :
126 : // iterate through the faces. If a face does not contain refEdge,
127 : // then we'll create a tetrahedron or a pyramid from the face to the
128 : // new vertex
129 0 : for(int i = 0; i < NUM_FACES; ++i){
130 0 : if(!FACE_CONTAINS_EDGE[i][refEdge]){
131 : // if the face is a triangle, then create a tetrahedron.
132 : // if it is a quadrilateral, then create a pyramid.
133 0 : const int* f = FACE_VRT_INDS[i];
134 :
135 0 : if(f[3] == -1){
136 : // create a tet
137 0 : inds[fi++] = GOID_TETRAHEDRON;
138 0 : inds[fi++] = f[0]; inds[fi++] = f[1];
139 0 : inds[fi++] = f[2]; inds[fi++] = nVrt;
140 : }
141 : else{
142 : // create a pyramid
143 0 : inds[fi++] = GOID_PYRAMID;
144 0 : inds[fi++] = f[0]; inds[fi++] = f[1];
145 0 : inds[fi++] = f[2]; inds[fi++] = f[3];
146 0 : inds[fi++] = nVrt;
147 : }
148 : }
149 : }
150 : }break;
151 :
152 0 : case 2: // multiple cases exist. Only some are directly supported.
153 : {
154 : // get the face in which the refEdges lie (if there is one at all)
155 0 : int refFaceInd = FACE_FROM_EDGES[refEdgeInds[0]][refEdgeInds[1]];
156 :
157 0 : if(refFaceInd != -1){
158 : // both refined edges lie in one face.
159 : // they either lie in the quad or in one of the triangles
160 0 : const int* f = FACE_VRT_INDS[refFaceInd];
161 0 : if(f[3] != -1){
162 : // a quad
163 : // we have to check whether the vertices lie at opposing
164 : // edges or at adjacent edges. We'll check corner-states for that.
165 : // if there is a corner with status 2, they lie on adjacent ones.
166 : int corner2 = -1;
167 : int corner2LocInd = -1;
168 0 : for(int i = 0; i < 4; ++i){
169 0 : if(cornerStatus[f[i]] == 2){
170 : corner2 = f[i];
171 : corner2LocInd = i;
172 : break;
173 : }
174 : }
175 :
176 0 : if(corner2 == -1){
177 : // opposing edges. Rotate the pyramid so that the edge
178 : // between v0 and v1 will not be splitted.
179 : int steps = 0;
180 : int firstEdgeInd = EDGE_FROM_VRTS[0][1];
181 0 : if(newEdgeVrts[firstEdgeInd] != 0)
182 : steps = 1;
183 :
184 : // the rotated pyramid
185 : int p[5];
186 0 : RotatePyramid(p, steps);
187 :
188 : // important vertex indices
189 0 : const int v1v2 = EDGE_FROM_VRTS[p[1]][p[2]] + E;
190 0 : const int v3v0 = EDGE_FROM_VRTS[p[3]][p[0]] + E;
191 :
192 : // now build two pyramids
193 : int& fi = fillCount;
194 : int* inds = newIndsOut;
195 :
196 0 : inds[fi++] = GOID_PYRAMID;
197 0 : inds[fi++] = p[0]; inds[fi++] = p[1];
198 0 : inds[fi++] = v1v2; inds[fi++] = v3v0;
199 0 : inds[fi++] = p[4];
200 :
201 0 : inds[fi++] = GOID_PYRAMID;
202 0 : inds[fi++] = v3v0; inds[fi++] = v1v2;
203 0 : inds[fi++] = p[2]; inds[fi++] = p[3];
204 0 : inds[fi++] = p[4];
205 : }
206 : else{
207 : // adjacent edges. Rotate the pyramid so that the corner is
208 : // at vertex 2
209 : int p[5];
210 0 : int steps = 2 - corner2LocInd;
211 0 : RotatePyramid(p, steps);
212 :
213 : // important vertex indices
214 0 : const int v1v2 = EDGE_FROM_VRTS[p[1]][p[2]] + E;
215 0 : const int v2v3 = EDGE_FROM_VRTS[p[2]][p[3]] + E;
216 :
217 : // now we have to create 4 tetrahedrons
218 : int& fi = fillCount;
219 : int* inds = newIndsOut;
220 :
221 0 : inds[fi++] = GOID_TETRAHEDRON;
222 0 : inds[fi++] = p[0]; inds[fi++] = p[1];
223 0 : inds[fi++] = v1v2; inds[fi++] = p[4];
224 :
225 0 : inds[fi++] = GOID_TETRAHEDRON;
226 0 : inds[fi++] = v1v2; inds[fi++] = p[2];
227 0 : inds[fi++] = v2v3; inds[fi++] = p[4];
228 :
229 0 : inds[fi++] = GOID_TETRAHEDRON;
230 0 : inds[fi++] = v2v3; inds[fi++] = p[3];
231 0 : inds[fi++] = p[0]; inds[fi++] = p[4];
232 :
233 0 : inds[fi++] = GOID_TETRAHEDRON;
234 0 : inds[fi++] = v1v2; inds[fi++] = v2v3;
235 0 : inds[fi++] = p[0]; inds[fi++] = p[4];
236 : }
237 : }
238 : else{
239 : // a tri.
240 : // only nice cuts are currently directly supported
241 0 : if(cornerStatus[4] == 2){
242 : // we'll create a prism and a pyramid. However, first
243 : // we'll rotate the pyramid, so that the cut lies in face 2
244 : int p[5];
245 0 : RotatePyramid(p, (2 - refFaceInd));
246 :
247 0 : int e1 = EDGE_FROM_VRTS[p[1]][p[4]] + E;
248 0 : int e2 = EDGE_FROM_VRTS[p[2]][p[4]] + E;
249 :
250 : int& fi = fillCount;
251 : int* inds = newIndsOut;
252 :
253 0 : inds[fi++] = GOID_PRISM;
254 0 : inds[fi++] = p[1]; inds[fi++] = p[0]; inds[fi++] = e1;
255 0 : inds[fi++] = p[2]; inds[fi++] = p[3]; inds[fi++] = e2;
256 :
257 0 : inds[fi++] = GOID_PYRAMID;
258 0 : inds[fi++] = p[0]; inds[fi++] = e1; inds[fi++] = e2;
259 0 : inds[fi++] = p[3]; inds[fi++] = p[4];
260 : }
261 : }
262 : }
263 : }break;
264 :
265 0 : case 3:
266 : {
267 : // currently we only directly support refinement of the base
268 : // get the face in which the refEdges lie (if there is one at all)
269 0 : int refFaceInd = FACE_FROM_EDGES[refEdgeInds[0]][refEdgeInds[1]];
270 0 : int tFaceInd = FACE_FROM_EDGES[refEdgeInds[1]][refEdgeInds[2]];
271 :
272 0 : if(refFaceInd == tFaceInd && refFaceInd != -1){
273 0 : const int* f = FACE_VRT_INDS[refFaceInd];
274 : // make sure that the face is a quadrilateral
275 0 : if(f[3] != -1){
276 : // rotate the pyramid, so that the edge v0v1 is not refined
277 : int freeEdge = -1;
278 0 : for(int i = 0; i < 4; ++i){
279 0 : int e = FACE_EDGE_INDS[refFaceInd][i];
280 : assert(e != -1);
281 0 : if(!newEdgeVrts[e]){
282 : freeEdge = e;
283 : break;
284 : }
285 : }
286 :
287 : assert(freeEdge != -1);
288 :
289 : // the rotated pyramid
290 : int p[5];
291 0 : RotatePyramid(p, -freeEdge);
292 :
293 : // important vertex indices
294 0 : const int v1v2 = EDGE_FROM_VRTS[p[1]][p[2]] + E;
295 0 : const int v2v3 = EDGE_FROM_VRTS[p[2]][p[3]] + E;
296 0 : const int v3v0 = EDGE_FROM_VRTS[p[3]][p[0]] + E;
297 :
298 : // now we have to create 1 pyramid and 3 tetrahedrons
299 : int& fi = fillCount;
300 : int* inds = newIndsOut;
301 :
302 0 : inds[fi++] = GOID_PYRAMID;
303 0 : inds[fi++] = p[0]; inds[fi++] = p[1]; inds[fi++] = v1v2;
304 0 : inds[fi++] = v3v0; inds[fi++] = p[4];
305 :
306 0 : inds[fi++] = GOID_TETRAHEDRON;
307 0 : inds[fi++] = v1v2; inds[fi++] = p[2];
308 0 : inds[fi++] = v2v3; inds[fi++] = p[4];
309 :
310 0 : inds[fi++] = GOID_TETRAHEDRON;
311 0 : inds[fi++] = v2v3; inds[fi++] = p[3];
312 0 : inds[fi++] = v3v0; inds[fi++] = p[4];
313 :
314 0 : inds[fi++] = GOID_TETRAHEDRON;
315 0 : inds[fi++] = v1v2; inds[fi++] = v2v3;
316 0 : inds[fi++] = v3v0; inds[fi++] = p[4];
317 : }
318 : }
319 : }break;
320 :
321 0 : case 4:
322 : {
323 : // we only treat the two most important cases here: either all new
324 : // vertices lie in the base or none of them does.
325 0 : if(cornerStatus[4] == 0){
326 : // all lie in the base
327 : // we require an additional face vertex here
328 : const int nVrt = FACE_FROM_VRTS[0][1][2] + F;
329 : const int v0v1 = EDGE_FROM_VRTS[0][1] + E;
330 : const int v1v2 = EDGE_FROM_VRTS[1][2] + E;
331 : const int v2v3 = EDGE_FROM_VRTS[2][3] + E;
332 : const int v3v0 = EDGE_FROM_VRTS[3][0] + E;
333 :
334 : // create 4 new pyramids
335 : int& fi = fillCount;
336 : int* inds = newIndsOut;
337 :
338 0 : inds[fi++] = GOID_PYRAMID;
339 0 : inds[fi++] = 0; inds[fi++] = v0v1; inds[fi++] = nVrt;
340 0 : inds[fi++] = v3v0; inds[fi++] = 4;
341 :
342 0 : inds[fi++] = GOID_PYRAMID;
343 0 : inds[fi++] = 1; inds[fi++] = v1v2; inds[fi++] = nVrt;
344 0 : inds[fi++] = v0v1; inds[fi++] = 4;
345 :
346 0 : inds[fi++] = GOID_PYRAMID;
347 0 : inds[fi++] = 2; inds[fi++] = v2v3; inds[fi++] = nVrt;
348 0 : inds[fi++] = v1v2; inds[fi++] = 4;
349 :
350 0 : inds[fi++] = GOID_PYRAMID;
351 0 : inds[fi++] = 3; inds[fi++] = v3v0; inds[fi++] = nVrt;
352 0 : inds[fi++] = v2v3; inds[fi++] = 4;
353 : }
354 0 : else if(cornerStatus[4] == 4){
355 : // all lie on the side edges
356 : const int e0 = EDGE_FROM_VRTS[0][4] + E;
357 : const int e1 = EDGE_FROM_VRTS[1][4] + E;
358 : const int e2 = EDGE_FROM_VRTS[2][4] + E;
359 : const int e3 = EDGE_FROM_VRTS[3][4] + E;
360 :
361 : // create a hexahedron and a new pyramid
362 : int& fi = fillCount;
363 : int* inds = newIndsOut;
364 :
365 0 : inds[fi++] = GOID_HEXAHEDRON;
366 0 : inds[fi++] = 0; inds[fi++] = 1;
367 0 : inds[fi++] = 2; inds[fi++] = 3;
368 0 : inds[fi++] = e0; inds[fi++] = e1;
369 0 : inds[fi++] = e2; inds[fi++] = e3;
370 :
371 0 : inds[fi++] = GOID_PYRAMID;
372 0 : inds[fi++] = e0; inds[fi++] = e1; inds[fi++] = e2;
373 0 : inds[fi++] = e3; inds[fi++] = 4;
374 : }
375 : }break;
376 :
377 0 : case 6:
378 : {
379 : // in the considered cases all 4 top edges and two opposing bottom edges
380 : // are marked
381 0 : if(cornerStatus[TOP_VERTEX] == 4){
382 : int rotSteps = -1;
383 0 : if(newEdgeVrts[BOTTOM_EDGE_INDS[0]] && newEdgeVrts[BOTTOM_EDGE_INDS[2]])
384 : rotSteps = 0;
385 0 : else if(newEdgeVrts[BOTTOM_EDGE_INDS[1]] && newEdgeVrts[BOTTOM_EDGE_INDS[3]])
386 : rotSteps = 1;
387 :
388 : if(rotSteps != -1){
389 : int vrts[NUM_VERTICES];
390 0 : RotatePyramid(vrts, rotSteps);
391 :
392 0 : const int e0 = EDGE_FROM_VRTS[vrts[0]][TOP_VERTEX] + E;
393 0 : const int e1 = EDGE_FROM_VRTS[vrts[1]][TOP_VERTEX] + E;
394 0 : const int e2 = EDGE_FROM_VRTS[vrts[2]][TOP_VERTEX] + E;
395 0 : const int e3 = EDGE_FROM_VRTS[vrts[3]][TOP_VERTEX] + E;
396 :
397 0 : const int v0v1 = EDGE_FROM_VRTS[vrts[0]][vrts[1]] + E;
398 0 : const int v2v3 = EDGE_FROM_VRTS[vrts[2]][vrts[3]] + E;
399 :
400 : // create 1 pyramid and 3 prisms
401 : int& fi = fillCount;
402 : int* inds = newIndsOut;
403 :
404 0 : inds[fi++] = GOID_PYRAMID;
405 0 : inds[fi++] = e0; inds[fi++] = e1; inds[fi++] = e2;
406 0 : inds[fi++] = e3; inds[fi++] = TOP_VERTEX;
407 :
408 0 : inds[fi++] = GOID_PRISM;
409 0 : inds[fi++] = vrts[0]; inds[fi++] = e0; inds[fi++] = v0v1;
410 0 : inds[fi++] = vrts[3]; inds[fi++] = e3; inds[fi++] = v2v3;
411 :
412 0 : inds[fi++] = GOID_PRISM;
413 0 : inds[fi++] = vrts[1]; inds[fi++] = v0v1; inds[fi++] = e1;
414 0 : inds[fi++] = vrts[2]; inds[fi++] = v2v3; inds[fi++] = e2;
415 :
416 0 : inds[fi++] = GOID_PRISM;
417 0 : inds[fi++] = v0v1; inds[fi++] = e0; inds[fi++] = e1;
418 0 : inds[fi++] = v2v3; inds[fi++] = e3; inds[fi++] = e2;
419 : }
420 : }
421 : }break;
422 :
423 0 : case 8: // regular refine
424 : {
425 : // we have to create 6 new pyramids and 4 tetrahedrons
426 : int& fi = fillCount;
427 : int* inds = newIndsOut;
428 0 : inds[fi++] = GOID_PYRAMID;
429 0 : inds[fi++] = 0; inds[fi++] = E;
430 0 : inds[fi++] = F; inds[fi++] = E + 3; inds[fi++] = E + 4;
431 :
432 0 : inds[fi++] = GOID_PYRAMID;
433 0 : inds[fi++] = E; inds[fi++] = 1;
434 0 : inds[fi++] = E + 1; inds[fi++] = F; inds[fi++] = E + 5;
435 :
436 0 : inds[fi++] = GOID_PYRAMID;
437 0 : inds[fi++] = F; inds[fi++] = E + 1;
438 0 : inds[fi++] = 2; inds[fi++] = E + 2; inds[fi++] = E + 6;
439 :
440 0 : inds[fi++] = GOID_PYRAMID;
441 0 : inds[fi++] = E + 3; inds[fi++] = F;
442 0 : inds[fi++] = E + 2; inds[fi++] = 3; inds[fi++] = E + 7;
443 :
444 0 : inds[fi++] = GOID_PYRAMID;
445 0 : inds[fi++] = E + 7; inds[fi++] = E + 6;
446 0 : inds[fi++] = E + 5; inds[fi++] = E + 4; inds[fi++] = F;
447 :
448 0 : inds[fi++] = GOID_PYRAMID;
449 0 : inds[fi++] = E + 4; inds[fi++] = E + 5;
450 0 : inds[fi++] = E + 6; inds[fi++] = E + 7; inds[fi++] = 4;
451 :
452 : // tetrahedrons
453 0 : inds[fi++] = GOID_TETRAHEDRON;
454 0 : inds[fi++] = E + 5; inds[fi++] = E;
455 0 : inds[fi++] = E + 4; inds[fi++] = F;
456 :
457 0 : inds[fi++] = GOID_TETRAHEDRON;
458 0 : inds[fi++] = E + 6; inds[fi++] = E + 1;
459 0 : inds[fi++] = E + 5; inds[fi++] = F;
460 :
461 0 : inds[fi++] = GOID_TETRAHEDRON;
462 0 : inds[fi++] = E + 7; inds[fi++] = E + 2;
463 0 : inds[fi++] = E + 6; inds[fi++] = F;
464 :
465 0 : inds[fi++] = GOID_TETRAHEDRON;
466 0 : inds[fi++] = E + 4; inds[fi++] = E + 3;
467 0 : inds[fi++] = E + 7; inds[fi++] = F;
468 :
469 0 : }break;
470 : }
471 :
472 :
473 0 : if(!snapPointsProcessed){
474 : UG_LOG("WARNING: Invalid or unsupported snap-point distribution detected. Ignoring snap-points for this element.\n");
475 : }
476 :
477 :
478 0 : if(fillCount == 0){
479 : // call recursive refine
480 0 : fillCount = shared_rules::RecursiveRefine(newIndsOut, newEdgeVrts,
481 : FACE_VRT_INDS, FACE_EDGE_INDS, NUM_VERTICES,
482 : NUM_EDGES, NUM_FACES);
483 :
484 : // the rule requires a new center vertex
485 0 : newCenterOut = true;
486 : }
487 :
488 0 : return fillCount;
489 : }
490 :
491 :
492 0 : bool IsRegularRefRule(const int edgeMarks)
493 : {
494 0 : return false;
495 : }
496 :
497 :
498 0 : int CollapseEdge (int* newIndsOut, int v0, int v1)
499 : {
500 : // get the edge which is connected by v0 and v1
501 0 : const int edgeInd = EDGE_FROM_VRTS[v0][v1];
502 :
503 : // if the edge doesn't exist or if one of the triangle-edges is collapsed,
504 : // no volume can be created
505 0 : if((edgeInd == -1) || (edgeInd > 3))
506 : return 0;
507 :
508 : // the index of the opposing edge of the base quadrilateral
509 0 : const int opEdgeInd = (edgeInd + 2) % 4;
510 :
511 : // the resulting volume is a tetrahedron
512 : int fi = 0;
513 0 : newIndsOut[fi++] = GOID_TETRAHEDRON;
514 0 : newIndsOut[fi++] = v0;
515 0 : newIndsOut[fi++] = EDGE_VRT_INDS[opEdgeInd][0];
516 0 : newIndsOut[fi++] = EDGE_VRT_INDS[opEdgeInd][1];
517 0 : newIndsOut[fi++] = TOP_VERTEX;
518 :
519 0 : return fi;
520 : }
521 :
522 : }// end of namespace
523 : }// end of namespace
|
