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 <algorithm>
35 : #include <iostream>
36 : #include "rule_util.h"
37 : #include "tetrahedron_rules.h"
38 : #include "octahedron_rules.h"
39 : #include "pyramid_rules.h"
40 : #include "hexahedron_rules.h"
41 :
42 : using namespace std;
43 :
44 : namespace ug{
45 : namespace shared_rules{
46 :
47 0 : int RecursiveRefine(int* newIndsOut, int* newEdgeVrts,
48 : const int faceVrtInds[][4],
49 : const int faceEdgeInds[][4],
50 : int numVrts, int numEdges, int numFaces)
51 : {
52 : // no refinement rule was found. We'll thus insert a new vertex in the
53 : // center and split the element into sub-elements.
54 : // we'll then recursively call this method for each sub-tetrahedron/-pyramid.
55 : // since all refinement rules for one side exist for those element types,
56 : // one recursion will always suffice.
57 :
58 : // a helper which avoids multiple recursion
59 : /*static bool dbgRecursionActive = false;
60 : assert(!dbgRecursionActive);
61 : dbgRecursionActive = true;*/
62 :
63 : // #define LIBGRID_RECURSIVE_REFINE__PRINT_DEBUG_INFO
64 : #ifdef LIBGRID_RECURSIVE_REFINE__PRINT_DEBUG_INFO
65 : static bool firstVisit = true;
66 : if(firstVisit){
67 : firstVisit = false;
68 : UG_LOG("DEBUGGING: LIBGRID_RECURSIVE_REFINE__PRINT_DEBUG_INFO enabled.\n");
69 : }
70 : UG_LOG("DEBUG: RecursiveRefine called on element with " << numVrts << " vertices.\n");
71 : UG_LOG(" Edge mark pattern: ");
72 : for(int i = 0; i < numEdges; ++i){
73 : cout << " " << newEdgeVrts[i] != NULL;
74 : }
75 : cout << endl;
76 : #endif
77 :
78 :
79 : int fillCount = 0;
80 :
81 : // we use constants that are chosen big enough to fit.
82 : int tmpInds[pyra_rules::MAX_NUM_INDS_OUT];
83 : int tmpNewEdgeVrts[12];
84 :
85 : // the indMap is used to map the generated temporary indices to the
86 : // indices of the main element.
87 0 : const int indMapSize = numVrts + numEdges + 1; // +1 for a quad-face
88 : int indMap[hex_rules::NUM_VERTICES + hex_rules::NUM_EDGES + 1];
89 : assert(indMapSize <= hex_rules::NUM_VERTICES + hex_rules::NUM_EDGES + 1);
90 :
91 0 : for(int i_face = 0; i_face < numFaces; ++i_face){
92 : // fill indMap and tmpNewEdgeVrts.
93 0 : const int* f = faceVrtInds[i_face];
94 : int numVrtsOfNewElem = 5;
95 0 : if(f[3] == -1)
96 : numVrtsOfNewElem = 4;
97 :
98 0 : for(int i_ind = 0; i_ind < indMapSize; ++i_ind)
99 0 : indMap[i_ind] = -1;
100 :
101 0 : for(int i = 0; i < 4; ++i){
102 0 : if(f[i] != -1)
103 0 : indMap[i] = f[i];
104 : }
105 :
106 0 : for(int i_edge = 0; i_edge < numEdges; ++i_edge)
107 0 : tmpNewEdgeVrts[i_edge] = 0;
108 :
109 0 : for(int i_edge = 0; i_edge < 4; ++i_edge){
110 0 : const int edgeInd = faceEdgeInds[i_face][i_edge];
111 0 : if(edgeInd != -1 && newEdgeVrts[edgeInd]){
112 0 : tmpNewEdgeVrts[i_edge] = 1;
113 0 : indMap[numVrtsOfNewElem + i_edge] = edgeInd + numVrts;
114 : }
115 : }
116 :
117 : // assign the index of the face vertex (only relevant for quad-faces) and
118 : // assign index of inner vertex to indMap and call refine.
119 0 : bool centerVrtRequired = false;
120 : int count = 0;
121 0 : if(f[3] == -1){
122 : //clog << "calling tet_rules::Refine\n";
123 0 : indMap[3] = numVrts + numEdges + numFaces;
124 0 : count = tet_rules::Refine(tmpInds, tmpNewEdgeVrts, centerVrtRequired);
125 : }
126 : else{
127 0 : indMap[4] = numVrts + numEdges + numFaces;
128 0 : indMap[pyra_rules::NUM_VERTICES + pyra_rules::NUM_EDGES] =
129 0 : numVrts + numEdges + i_face;
130 : //clog << "calling pyra_rules::Refine\n";
131 0 : count = pyra_rules::Refine(tmpInds, tmpNewEdgeVrts, centerVrtRequired);
132 : }
133 :
134 : // No center vertex should be created during this recudsion.
135 : assert(!centerVrtRequired);
136 : /*
137 : clog << "tmpInds:";
138 : for(int i = 0; i < count; ++i)
139 : clog << " " << tmpInds[i];
140 : clog << endl;
141 :
142 : clog << "indMap:";
143 : for(int i = 0; i < indMapSize; ++i)
144 : clog << " " << indMap[i];
145 : clog << endl;
146 : */
147 : // copy the new indices to the output array
148 0 : for(int i = 0; i < count;){
149 0 : int numElemInds = tmpInds[i];
150 0 : newIndsOut[fillCount + i] = numElemInds;
151 0 : ++i;
152 :
153 0 : for(int j = 0; j < numElemInds; ++j, ++i){
154 : assert(indMap[tmpInds[i]] != -1);
155 0 : newIndsOut[fillCount + i] = indMap[tmpInds[i]];
156 : }
157 : }
158 :
159 0 : fillCount += count;
160 : }
161 :
162 : //dbgRecursionActive = false;
163 0 : return fillCount;
164 : }
165 :
166 : }// end of namespace
167 : }// end of namespace
|
