1 /*--------------------------------------------------------------------------*/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /* file: traverse_nr_2d.c */
7 /* */
8 /* description: */
9 /* nonrecursive mesh traversal, 2d routines */
10 /* */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Alfred Schmidt */
14 /* Zentrum fuer Technomathematik */
15 /* Fachbereich 3 Mathematik/Informatik */
16 /* Universitaet Bremen */
17 /* Bibliothekstr. 2 */
18 /* D-28359 Bremen, Germany */
19 /* */
20 /* Kunibert G. Siebert */
21 /* Institut fuer Mathematik */
22 /* Universitaet Augsburg */
23 /* Universitaetsstr. 14 */
24 /* D-86159 Augsburg, Germany */
25 /* */
26 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
27 /* */
28 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
29 /* */
30 /*--------------------------------------------------------------------------*/
31
32 static int coarse_nb_2d[3][3] = {{-2,-2,-2}, {2,-1,1}, {-1,2,0}};
33 /* father.neigh[coarse_nb_2d[i][j]] == child[i-1].neigh[j] */
34
35 /*--------------------------------------------------------------------------*/
36 /* traverse_neighbour_2d: */
37 /* ------------------- */
38 /* walk through hierarchy tree and look for a neighbour */
39 /*--------------------------------------------------------------------------*/
40
41
traverse_neighbour_2d(TRAVERSE_STACK * stack,const EL_INFO * elinfo_old,int neighbour)42 static const EL_INFO *traverse_neighbour_2d(TRAVERSE_STACK *stack,
43 const EL_INFO *elinfo_old,
44 int neighbour)
45 {
46 FUNCNAME("traverse_neighbour_2d");
47 EL *el, *sav_el;
48 EL_INFO *old_elinfo, *elinfo;
49 int i, nb, opp_vertex, stack2_used;
50 int sav_index, sav_neighbour = neighbour;
51 int parity = 0;
52
53 DEBUG_TEST_EXIT(stack->stack_used > 0, "no current element");
54 DEBUG_TEST_EXIT(stack->traverse_flags & CALL_LEAF_EL,
55 "invalid traverse_fill_flag=%d", stack->traverse_flags);
56 DEBUG_TEST_EXIT(elinfo_old == stack->elinfo_stack+stack->stack_used,
57 "invalid old elinfo");
58 DEBUG_TEST_FLAG(FILL_NEIGH, stack->elinfo_stack+stack->stack_used);
59
60 el = stack->elinfo_stack[stack->stack_used].el;
61 sav_index = INDEX(el);
62 sav_el = el;
63
64 /* first, goto to leaf level, if necessary... */
65 if ((el->child[0]) && (neighbour < 2)) {
66 if (stack->stack_used >= stack->stack_size-1)
67 __AI_enlarge_traverse_stack(stack);
68 i = 1 - neighbour;
69 fill_elinfo(i, stack->fill_flag, stack->elinfo_stack+stack->stack_used,
70 stack->elinfo_stack+stack->stack_used+1);
71 stack->info_stack[stack->stack_used] = i + 1;
72 stack->stack_used++;
73 neighbour = 2;
74 }
75
76 /* save information about current element and its position in the tree */
77 stack->save_traverse_mel = stack->traverse_mel;
78 stack->save_stack_used = stack->stack_used;
79 for (i=0; i<=stack->stack_used; i++)
80 stack->save_info_stack[i] = stack->info_stack[i];
81 for (i=0; i<=stack->stack_used; i++)
82 stack->save_elinfo_stack[i] = stack->elinfo_stack[i];
83 old_elinfo = stack->save_elinfo_stack+stack->stack_used;
84 opp_vertex = old_elinfo->opp_vertex[neighbour];
85
86 /*--------------------------------------------------------------------------*/
87 /* First phase: go up in tree until we can go down again. */
88 /* */
89 /* During this first phase, nb is the neighbour index which points from an */
90 /* element of the OLD hierarchy branch to the NEW branch */
91 /*--------------------------------------------------------------------------*/
92
93 nb = neighbour;
94
95 while (stack->stack_used > 1)
96 {
97 stack->stack_used--;
98 nb = coarse_nb_2d[stack->info_stack[stack->stack_used]][nb];
99 if (nb == -1) break;
100 DEBUG_TEST_EXIT(nb >= 0, "invalid coarse_nb_2d %d\n",nb);
101 }
102
103 /*--------------------------------------------------------------------------*/
104 /* Now, goto neighbouring element at the local hierarchy entry */
105 /* This is either a macro element neighbour or the other child of parent. */
106 /* initialize nb for second phase (see below) */
107 /*--------------------------------------------------------------------------*/
108
109 if (nb >= 0) { /* go to macro element neighbour */
110 const MACRO_EL *mel_old;
111
112 if (stack->traverse_mel->neigh[nb] == NULL) {
113 return NULL;
114 }
115
116 if ((nb < 2) && (stack->save_stack_used > 1)) {
117 stack2_used = 2; /* go down one level in OLD hierarchy */
118 } else {
119 stack2_used = 1;
120 }
121
122 mel_old = stack->traverse_mel;
123 i = mel_old->opp_vertex[nb];
124 stack->traverse_mel = stack->traverse_mel->neigh[nb];
125
126 if (mel_old->neigh_vertices[nb][0] != -1) {
127 if (mel_old->neigh_vertices[nb][0] == ((i+1) % N_VERTICES_2D)) {
128 parity = 1;
129 }
130 } else if (mel_old->coord[(nb + 1) % N_VERTICES_2D]
131 ==
132 stack->traverse_mel->coord[(i+1) % N_VERTICES_2D]) {
133 /* neighbour has different orientation */
134 parity = 1;
135 }
136
137 nb = i;
138
139 stack->stack_used = 1;
140 fill_macro_info(stack->traverse_mesh, stack->traverse_mel,
141 stack->elinfo_stack+stack->stack_used);
142 stack->info_stack[stack->stack_used] = 0;
143
144 } else { /* goto other child */
145
146 stack2_used = stack->stack_used + 1;
147 if (stack->save_stack_used > stack2_used) {
148 stack2_used++; /* go down one level in OLD hierarchy */
149 }
150 /* elinfo2 = stack->save_elinfo_stack+stack2_used; */
151 /* el2 = elinfo2->el; */
152
153
154 if (stack->stack_used >= stack->stack_size-1)
155 __AI_enlarge_traverse_stack(stack);
156 i = 2 - stack->info_stack[stack->stack_used];
157 stack->info_stack[stack->stack_used] = i+1;
158 fill_elinfo(i, stack->fill_flag, stack->elinfo_stack+stack->stack_used,
159 stack->elinfo_stack+stack->stack_used+1);
160 stack->stack_used++;
161 nb = 1-i;
162 }
163
164 /*--------------------------------------------------------------------------*/
165 /* Second phase: go down in a new hierarchy branch until leaf level. */
166 /* Now, nb is the neighbour index which points from an element of the */
167 /* NEW hierarchy branch to the OLD branch. */
168 /*--------------------------------------------------------------------------*/
169
170 elinfo = stack->elinfo_stack+stack->stack_used;
171 el = elinfo->el;
172
173 while(el->child[0]) {
174
175 if (nb < 2) { /* go down one level in hierarchy */
176 if (stack->stack_used >= stack->stack_size-1)
177 __AI_enlarge_traverse_stack(stack);
178 fill_elinfo(1-nb, stack->fill_flag,
179 stack->elinfo_stack+stack->stack_used,
180 stack->elinfo_stack+stack->stack_used+1);
181 stack->info_stack[stack->stack_used] = 2-nb;
182 stack->stack_used++;
183 nb = 2;
184 }
185
186 if (stack->save_stack_used > stack2_used) { /* `refine' both el and el2 */
187 DEBUG_TEST_EXIT(el->child[0], "invalid NEW refinement?");
188
189 /* use child i, neighbour of el2->child[nb-1] */
190 i = 2 - stack->save_info_stack[stack2_used];
191 DEBUG_TEST_EXIT(i < 2, "invalid OLD refinement?");
192 if (parity)
193 i = 1 - i;
194 stack->info_stack[stack->stack_used] = i+1;
195 fill_elinfo(i, stack->fill_flag,
196 stack->elinfo_stack+stack->stack_used,
197 stack->elinfo_stack+stack->stack_used+1);
198 stack->stack_used++;
199 nb = i;
200
201 elinfo = stack->elinfo_stack+stack->stack_used;
202 el = elinfo->el;
203
204 stack2_used++;
205 if (stack->save_stack_used > stack2_used) {
206 stack2_used++; /* go down one level in OLD hierarchy */
207 }
208 /* elinfo2 = stack->save_elinfo_stack+stack2_used; */
209 /* el2 = elinfo2->el; */
210
211 } else { /* now we're done... */
212 elinfo = stack->elinfo_stack+stack->stack_used;
213 el = elinfo->el;
214 }
215 }
216
217
218 if (elinfo->neigh[opp_vertex] != old_elinfo->el) {
219 MSG(" looking for neighbour %d of element %d at %p\n",
220 neighbour, INDEX(old_elinfo->el), old_elinfo->el);
221 MSG(" originally: neighbour %d of element %d at %p\n",
222 sav_neighbour, sav_index, sav_el);
223 MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
224 INDEX(elinfo->el), elinfo->el, opp_vertex,
225 INDEX(elinfo->neigh[opp_vertex]));
226 DEBUG_TEST_EXIT(elinfo->neigh[opp_vertex] == old_elinfo->el,
227 "didn't succeed !?!?!?");
228 }
229 if (elinfo->el->child[0])
230 {
231 MSG(" looking for neighbour %d of element %d at %p\n",
232 neighbour, INDEX(old_elinfo->el), old_elinfo->el);
233 MSG(" originally: neighbour %d of element %d at %p\n",
234 sav_neighbour, sav_index, sav_el);
235 MSG(" got element %d at %p with opp_vertex %d neigh %d\n",
236 INDEX(elinfo->el), elinfo->el, opp_vertex,
237 INDEX(elinfo->neigh[opp_vertex]));
238 MSG("got no leaf element\n");
239 WAIT_REALLY;
240 }
241
242 elinfo->el_geom_cache.current_el = NULL;
243
244 return elinfo;
245 }
246
247