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 "rule_util.h"
35 : #include "hexahedron_rules.h"
36 : #include "grid_object_ids.h"
37 :
38 : namespace ug{
39 : namespace hex_rules{
40 :
41 : /// Output are the vertices of a rotated quadrilateral.
42 0 : void RotateQuad(int vrtsOut[4], const int quad[4], 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] = quad[i];
49 0 : }
50 :
51 0 : int Refine(int* newIndsOut, int* newEdgeVrts, bool& newCenterOut, vector3*,
52 : bool* isSnapPoint)
53 : {
54 0 : newCenterOut = false;
55 : // If a refinement rule is not implemented, fillCount will stay at 0.
56 : // Before returning, we'll check, whether fillCount is at 0 and will
57 : // perform refinement with an additional vertex in this case.
58 :
59 : // corner status is used to mark vertices, which are connected to refined edges
60 : // the value tells, how many associated edges are refined.
61 0 : int cornerStatus[8] = {0, 0, 0, 0, 0, 0, 0, 0};
62 :
63 : // here we'll store the index of each edge, which will be refined.
64 : // Size of the valid sub-array is numNewVrts, which is calculated below.
65 : int refEdgeInds[NUM_EDGES];
66 :
67 : // count the number of new vertices and fill newVrtEdgeInds
68 : int numNewVrts = 0;
69 0 : for(int i = 0; i < NUM_EDGES; ++i){
70 0 : if(newEdgeVrts[i]){
71 0 : refEdgeInds[numNewVrts] = i;
72 0 : ++numNewVrts;
73 :
74 : // adjust corner status of associated vertices
75 0 : const int* evi = EDGE_VRT_INDS[i];
76 0 : ++cornerStatus[evi[0]];
77 0 : ++cornerStatus[evi[1]];
78 : }
79 : }
80 :
81 : // snap-point handling
82 : int numSnapPoints = 0;
83 : // int snapPointIndex[NUM_VERTICES];
84 0 : if(isSnapPoint){
85 0 : for(int i = 0; i < NUM_VERTICES; ++i){
86 0 : if(isSnapPoint[i]){
87 : // snapPointIndex[numSnapPoints] = i;
88 0 : ++numSnapPoints;
89 : }
90 : }
91 : }
92 :
93 0 : bool snapPointsProcessed = (numSnapPoints == 0);
94 :
95 : // the fillCount tells how much data has already been written to newIndsOut.
96 : int fillCount = 0;
97 :
98 : // convenience - indices where new edge-vrts, new face-vrts and new vol-vrts begin.
99 : const int E = NUM_VERTICES;
100 : const int F = NUM_VERTICES + NUM_EDGES;
101 : const int V = NUM_VERTICES + NUM_EDGES + NUM_FACES;
102 :
103 : // depending on the number of new vertices, we will now apply different
104 : // refinement rules. Further distinction may have to be done.
105 0 : switch(numNewVrts){
106 0 : case 0:
107 : {
108 : // simply put the default hexahedron back to newIndsOut
109 0 : newIndsOut[fillCount++] = GOID_HEXAHEDRON;
110 0 : newIndsOut[fillCount++] = 0;
111 0 : newIndsOut[fillCount++] = 1;
112 0 : newIndsOut[fillCount++] = 2;
113 0 : newIndsOut[fillCount++] = 3;
114 0 : newIndsOut[fillCount++] = 4;
115 0 : newIndsOut[fillCount++] = 5;
116 0 : newIndsOut[fillCount++] = 6;
117 0 : newIndsOut[fillCount++] = 7;
118 0 : }break;
119 :
120 0 : case 1:
121 : {
122 : // we'll iterate over all faces. If a face does not contain the
123 : // refined edge, then we'll create a pyramid from it.
124 0 : const int refEdge = refEdgeInds[0];
125 0 : const int nVrt = refEdge + E;
126 :
127 : int& fi = fillCount;
128 : int* inds = newIndsOut;
129 :
130 0 : for(int i = 0; i < NUM_FACES; ++i){
131 0 : if(!FACE_CONTAINS_EDGE[i][refEdge]){
132 0 : const int* f = FACE_VRT_INDS[i];
133 0 : inds[fi++] = GOID_PYRAMID;
134 0 : inds[fi++] = f[0]; inds[fi++] = f[1];
135 0 : inds[fi++] = f[2]; inds[fi++] = f[3];
136 0 : inds[fi++] = nVrt;
137 : }
138 : }
139 : }break;
140 :
141 0 : case 2:
142 : {
143 : // currently only the case where both new vertices lie on opposed
144 : // edges of one side is supported.
145 : // We thus try to get that face.
146 0 : int refFaceInd = FACE_FROM_EDGES[refEdgeInds[0]][refEdgeInds[1]];
147 0 : if(refFaceInd == -1)
148 : break;
149 :
150 : // make sure that corners of refEdge0 both have corner-state 1
151 0 : const int* refEdge0 = EDGE_VRT_INDS[refEdgeInds[0]];
152 0 : if(cornerStatus[refEdge0[0]] != 1 || cornerStatus[refEdge0[1]] != 1)
153 : break;
154 :
155 : // find the first face != refFaceInd, which contains refEdge0
156 : int faceInd = -1;
157 0 : for(int i = 0; i < NUM_FACES; ++i){
158 0 : if(FACE_CONTAINS_EDGE[i][refEdgeInds[0]] && i != refFaceInd){
159 : faceInd = i;
160 : break;
161 : }
162 : }
163 :
164 : assert(faceInd != -1);
165 :
166 : // the face and its opposing face
167 0 : const int* tf = FACE_VRT_INDS[faceInd];
168 0 : const int* tof = OPPOSED_FACE_VRT_INDS[faceInd];
169 :
170 : // we have to rotate the faces, so that the rules can be easily applied
171 : int localEdgeInd = -1;
172 0 : for(int i = 0; i < 4; ++i){
173 0 : if(newEdgeVrts[EDGE_FROM_VRTS[tf[i]][tf[(i + 1) % 4]]]){
174 : localEdgeInd = i;
175 : break;
176 : }
177 : }
178 :
179 : assert(localEdgeInd != -1);
180 :
181 : // the local refEdge shall lie between vrts 1, 2
182 : int f[4];
183 : int of[4];
184 0 : RotateQuad(f, tf, 1 - localEdgeInd);
185 0 : RotateQuad(of, tof, 1 - localEdgeInd);
186 :
187 : // create three new prisms
188 0 : const int v0 = EDGE_FROM_VRTS[f[1]][f[2]] + E;
189 0 : const int v1 = EDGE_FROM_VRTS[of[1]][of[2]] + E;
190 :
191 : int& fi = fillCount;
192 : int* inds = newIndsOut;
193 :
194 0 : if(numSnapPoints == 2){
195 : // two cases are supported.
196 0 : if(isSnapPoint[f[0]] && isSnapPoint[of[0]]){
197 0 : inds[fi++] = GOID_PRISM;
198 0 : inds[fi++] = f[0]; inds[fi++] = f[1]; inds[fi++] = v0;
199 0 : inds[fi++] = of[0]; inds[fi++] = of[1]; inds[fi++] = v1;
200 :
201 0 : inds[fi++] = GOID_HEXAHEDRON;
202 0 : inds[fi++] = f[0]; inds[fi++] = v0; inds[fi++] = f[2]; inds[fi++] = f[3];
203 0 : inds[fi++] = of[0]; inds[fi++] = v1; inds[fi++] = of[2]; inds[fi++] = of[3];
204 0 : snapPointsProcessed = true;
205 : }
206 0 : else if(isSnapPoint[f[3]] && isSnapPoint[of[3]]){
207 0 : inds[fi++] = GOID_PRISM;
208 0 : inds[fi++] = f[2]; inds[fi++] = f[3]; inds[fi++] = v0;
209 0 : inds[fi++] = of[2]; inds[fi++] = of[3]; inds[fi++] = v1;
210 :
211 0 : inds[fi++] = GOID_HEXAHEDRON;
212 0 : inds[fi++] = f[0]; inds[fi++] = f[1]; inds[fi++] = v0; inds[fi++] = f[3];
213 0 : inds[fi++] = of[0]; inds[fi++] = of[1]; inds[fi++] = v1; inds[fi++] = of[3];
214 : snapPointsProcessed = true;
215 : }
216 : }
217 :
218 0 : if(numSnapPoints == 0 || !snapPointsProcessed){
219 0 : inds[fi++] = GOID_PRISM;
220 0 : inds[fi++] = f[0]; inds[fi++] = f[1]; inds[fi++] = v0;
221 0 : inds[fi++] = of[0]; inds[fi++] = of[1]; inds[fi++] = v1;
222 :
223 0 : inds[fi++] = GOID_PRISM;
224 0 : inds[fi++] = f[0]; inds[fi++] = v0; inds[fi++] = f[3];
225 0 : inds[fi++] = of[0]; inds[fi++] = v1; inds[fi++] = of[3];
226 :
227 0 : inds[fi++] = GOID_PRISM;
228 0 : inds[fi++] = f[2]; inds[fi++] = f[3]; inds[fi++] = v0;
229 0 : inds[fi++] = of[2]; inds[fi++] = of[3]; inds[fi++] = v1;
230 : }
231 0 : }break;
232 :
233 : case 4:
234 : {
235 : // we only support nice regular cuts with two resulting hexahedrons.
236 : // get the face which contains refEdge0
237 : int faceInd = -1;
238 0 : for(int i = 0; i < NUM_FACES; ++i){
239 0 : if(FACE_CONTAINS_EDGE[i][refEdgeInds[0]]){
240 : faceInd = i;
241 : break;
242 : }
243 : }
244 : assert(faceInd != -1);
245 :
246 0 : const int* tf = FACE_VRT_INDS[faceInd];
247 0 : const int* tof = OPPOSED_FACE_VRT_INDS[faceInd];
248 :
249 : // get the local index of refEdgeInds[0]
250 : int localEdgeInd = -1;
251 0 : for(int i = 0; i < 4; ++i){
252 0 : if(newEdgeVrts[EDGE_FROM_VRTS[tf[i]][tf[(i + 1) % 4]]]){
253 : localEdgeInd = i;
254 : break;
255 : }
256 : }
257 : assert(localEdgeInd != -1);
258 :
259 : // Check whether the opposed edge in tf will be refined, too
260 0 : if(newEdgeVrts[EDGE_FROM_VRTS[tf[(localEdgeInd + 2) % 4]]
261 0 : [tf[(localEdgeInd + 3) % 4]]])
262 : {
263 : // make sure that those are the same local indices as in tof
264 0 : if(!newEdgeVrts[EDGE_FROM_VRTS[tof[(localEdgeInd)]]
265 0 : [tof[(localEdgeInd + 1) % 4]]]
266 0 : || !newEdgeVrts[EDGE_FROM_VRTS[tof[(localEdgeInd + 2) % 4]]
267 0 : [tof[(localEdgeInd + 3) % 4]]])
268 : {
269 : // they are not. we can't handle that case.
270 : break;
271 : }
272 : }
273 : else
274 : break;
275 :
276 : // we rotate the faces, so that the refined edges always lie
277 : // between vertices 1,2 and 3,0
278 : int f[4];
279 : int of[4];
280 :
281 0 : RotateQuad(f, tf, 1 - localEdgeInd);
282 0 : RotateQuad(of, tof, 1 - localEdgeInd);
283 :
284 : int v0, v1, ov0, ov1;
285 :
286 : // get the new vertex indices of new vertices
287 0 : v0 = EDGE_FROM_VRTS[f[1]][f[2]] + E;
288 0 : v1 = EDGE_FROM_VRTS[f[3]][f[0]] + E;
289 0 : ov0 = EDGE_FROM_VRTS[of[1]][of[2]] + E;
290 0 : ov1 = EDGE_FROM_VRTS[of[3]][of[0]] + E;
291 :
292 : // create 2 new hexahedrons
293 : int& fi = fillCount;
294 : int* inds = newIndsOut;
295 :
296 0 : inds[fi++] = GOID_HEXAHEDRON;
297 0 : inds[fi++] = f[0]; inds[fi++] = f[1];
298 0 : inds[fi++] = v0; inds[fi++] = v1;
299 0 : inds[fi++] = of[0]; inds[fi++] = of[1];
300 0 : inds[fi++] = ov0; inds[fi++] = ov1;
301 :
302 0 : inds[fi++] = GOID_HEXAHEDRON;
303 0 : inds[fi++] = f[2]; inds[fi++] = f[3];
304 0 : inds[fi++] = v1; inds[fi++] = v0;
305 0 : inds[fi++] = of[2]; inds[fi++] = of[3];
306 0 : inds[fi++] = ov1; inds[fi++] = ov0;
307 0 : }break;
308 :
309 : case 8:
310 : {
311 : // we currently only support the case where two opposing faces
312 : // are completly refined.
313 : // iterate through all faces, until two are found whose edges are all
314 : // to be refined.
315 : int refFaceInd[2] = {-1, -1};
316 0 : for(int i = 0; i < NUM_FACES; ++i){
317 0 : const int* e = FACE_EDGE_INDS[i];
318 : int j = 0;
319 0 : for(;j < 4; ++j){
320 0 : if(!newEdgeVrts[e[j]])
321 : break;
322 : }
323 :
324 0 : if(j == 4){
325 0 : if(refFaceInd[0] == -1)
326 : refFaceInd[0] = i;
327 : else{
328 : assert(refFaceInd[1] == -1);
329 : refFaceInd[1] = i;
330 : }
331 : }
332 : }
333 :
334 : // if two faces were found and if they are opposed faces, then
335 : // we can create the four new hexahedrons
336 0 : if(refFaceInd[0] != -1 && refFaceInd[1] != -1
337 0 : && refFaceInd[1] == OPPOSED_FACE[refFaceInd[0]])
338 : {
339 0 : const int* f = FACE_VRT_INDS[refFaceInd[0]];
340 0 : const int* of = OPPOSED_FACE_VRT_INDS[refFaceInd[0]];
341 :
342 0 : const int e0 = EDGE_FROM_VRTS[f[0]][f[1]] + E;
343 0 : const int e1 = EDGE_FROM_VRTS[f[1]][f[2]] + E;
344 0 : const int e2 = EDGE_FROM_VRTS[f[2]][f[3]] + E;
345 0 : const int e3 = EDGE_FROM_VRTS[f[3]][f[0]] + E;
346 :
347 0 : const int oe0 = EDGE_FROM_VRTS[of[0]][of[1]] + E;
348 0 : const int oe1 = EDGE_FROM_VRTS[of[1]][of[2]] + E;
349 0 : const int oe2 = EDGE_FROM_VRTS[of[2]][of[3]] + E;
350 0 : const int oe3 = EDGE_FROM_VRTS[of[3]][of[0]] + E;
351 :
352 0 : const int fvrt = refFaceInd[0] + F;
353 0 : const int ofvrt = refFaceInd[1] + F;
354 :
355 : // create 4 new hexahedrons
356 : int& fi = fillCount;
357 : int* inds = newIndsOut;
358 :
359 0 : inds[fi++] = GOID_HEXAHEDRON;
360 0 : inds[fi++] = f[0]; inds[fi++] = e0;
361 0 : inds[fi++] = fvrt; inds[fi++] = e3;
362 0 : inds[fi++] = of[0]; inds[fi++] = oe0;
363 0 : inds[fi++] = ofvrt; inds[fi++] = oe3;
364 :
365 0 : inds[fi++] = GOID_HEXAHEDRON;
366 0 : inds[fi++] = f[1]; inds[fi++] = e1;
367 0 : inds[fi++] = fvrt; inds[fi++] = e0;
368 0 : inds[fi++] = of[1]; inds[fi++] = oe1;
369 0 : inds[fi++] = ofvrt; inds[fi++] = oe0;
370 :
371 0 : inds[fi++] = GOID_HEXAHEDRON;
372 0 : inds[fi++] = f[2]; inds[fi++] = e2;
373 0 : inds[fi++] = fvrt; inds[fi++] = e1;
374 0 : inds[fi++] = of[2]; inds[fi++] = oe2;
375 0 : inds[fi++] = ofvrt; inds[fi++] = oe1;
376 :
377 0 : inds[fi++] = GOID_HEXAHEDRON;
378 0 : inds[fi++] = f[3]; inds[fi++] = e3;
379 0 : inds[fi++] = fvrt; inds[fi++] = e2;
380 0 : inds[fi++] = of[3]; inds[fi++] = oe3;
381 0 : inds[fi++] = ofvrt; inds[fi++] = oe2;
382 : }
383 0 : }break;
384 :
385 0 : case 12: // regular refine
386 : {
387 : // we have to create 8 new hexahedrons
388 : int& fi = fillCount;
389 : int* inds = newIndsOut;
390 0 : inds[fi++] = GOID_HEXAHEDRON;
391 0 : inds[fi++] = 0; inds[fi++] = E;
392 0 : inds[fi++] = F; inds[fi++] = E + 3;
393 0 : inds[fi++] = E + 4; inds[fi++] = F + 1;
394 0 : inds[fi++] = V; inds[fi++] = F + 4;
395 :
396 0 : inds[fi++] = GOID_HEXAHEDRON;
397 0 : inds[fi++] = E; inds[fi++] = 1;
398 0 : inds[fi++] = E + 1; inds[fi++] = F;
399 0 : inds[fi++] = F + 1; inds[fi++] = E + 5;
400 0 : inds[fi++] = F + 2; inds[fi++] = V;
401 :
402 0 : inds[fi++] = GOID_HEXAHEDRON;
403 0 : inds[fi++] = F; inds[fi++] = E + 1;
404 0 : inds[fi++] = 2; inds[fi++] = E + 2;
405 0 : inds[fi++] = V; inds[fi++] = F + 2;
406 0 : inds[fi++] = E + 6; inds[fi++] = F + 3;
407 :
408 0 : inds[fi++] = GOID_HEXAHEDRON;
409 0 : inds[fi++] = E + 3; inds[fi++] = F;
410 0 : inds[fi++] = E + 2; inds[fi++] = 3;
411 0 : inds[fi++] = F + 4; inds[fi++] = V;
412 0 : inds[fi++] = F + 3; inds[fi++] = E + 7;
413 :
414 0 : inds[fi++] = GOID_HEXAHEDRON;
415 0 : inds[fi++] = E + 4; inds[fi++] = F + 1;
416 0 : inds[fi++] = V; inds[fi++] = F + 4;
417 0 : inds[fi++] = 4; inds[fi++] = E + 8;
418 0 : inds[fi++] = F + 5; inds[fi++] = E + 11;
419 :
420 0 : inds[fi++] = GOID_HEXAHEDRON;
421 0 : inds[fi++] = F + 1; inds[fi++] = E + 5;
422 0 : inds[fi++] = F + 2; inds[fi++] = V;
423 0 : inds[fi++] = E + 8; inds[fi++] = 5;
424 0 : inds[fi++] = E + 9; inds[fi++] = F + 5;
425 :
426 0 : inds[fi++] = GOID_HEXAHEDRON;
427 0 : inds[fi++] = V; inds[fi++] = F + 2;
428 0 : inds[fi++] = E + 6; inds[fi++] = F + 3;
429 0 : inds[fi++] = F + 5; inds[fi++] = E + 9;
430 0 : inds[fi++] = 6; inds[fi++] = E + 10;
431 :
432 0 : inds[fi++] = GOID_HEXAHEDRON;
433 0 : inds[fi++] = F + 4; inds[fi++] = V;
434 0 : inds[fi++] = F + 3; inds[fi++] = E + 7;
435 0 : inds[fi++] = E + 11; inds[fi++] = F + 5;
436 0 : inds[fi++] = E + 10; inds[fi++] = 7;
437 :
438 : // the rule requires a new center vertex
439 0 : newCenterOut = true;
440 0 : }break;
441 : }
442 :
443 0 : if(!snapPointsProcessed){
444 : UG_LOG("WARNING: Invalid or unsupported snap-point distribution detected. Ignoring snap-points for this element.\n");
445 : }
446 :
447 0 : if(fillCount == 0){
448 : // call recursive refine
449 0 : fillCount = shared_rules::RecursiveRefine(newIndsOut, newEdgeVrts,
450 : FACE_VRT_INDS, FACE_EDGE_INDS, NUM_VERTICES,
451 : NUM_EDGES, NUM_FACES);
452 :
453 : // the rule requires a new center vertex
454 0 : newCenterOut = true;
455 : }
456 :
457 0 : return fillCount;
458 : }
459 :
460 0 : bool IsRegularRefRule(const int edgeMarks)
461 : {
462 : // static const int edges[3][4] = {{0, 2, 8, 10}, {4, 5, 6, 7}, {1, 3, 9, 11}};
463 : static const int allEdges[3] = { 1285, //010100000101}
464 : 240, //000011110000
465 : 2570}; //101000001010
466 : int clearedMarks = 0;
467 0 : for(int i = 0; i < 3; ++i){
468 0 : int t = edgeMarks & allEdges[i];
469 0 : if(t != 0 && t != allEdges[i])
470 : return false;
471 0 : clearedMarks |= t;
472 : }
473 :
474 0 : return clearedMarks != 0;
475 : }
476 :
477 : }// end of namespace
478 : }// end of namespace
|
