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: macro.c */
7 /* */
8 /* */
9 /* description: dimension dependent part of reading/writing macro */
10 /* triangulations for 3d */
11 /* */
12 /*--------------------------------------------------------------------------*/
13 /* */
14 /* authors: Alfred Schmidt */
15 /* Zentrum fuer Technomathematik */
16 /* Fachbereich 3 Mathematik/Informatik */
17 /* Universitaet Bremen */
18 /* Bibliothekstr. 2 */
19 /* D-28359 Bremen, Germany */
20 /* */
21 /* Kunibert G. Siebert */
22 /* Institut fuer Mathematik */
23 /* Universitaet Augsburg */
24 /* Universitaetsstr. 14 */
25 /* D-86159 Augsburg, Germany */
26 /* */
27 /* Daniel Koester */
28 /* Institut fuer Mathematik */
29 /* Albert-Ludwigs-Universitaet Freiburg */
30 /* Hermann-Herder-Str. 10 */
31 /* D-79104 Freiburg */
32 /* */
33 /* Claus-Justus Heine */
34 /* Abteilung fuer Angewandte Mathematik */
35 /* Albert-Ludwigs-Universitaet Freiburg */
36 /* Hermann-Herder-Str. 10 */
37 /* D-79104 Freiburg im Breisgau, Germany */
38 /* */
39 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
40 /* */
41 /* (c) by A. Schmidt and K.G. Siebert (1996-2003), */
42 /* C.-J. Heine (2006-2007) */
43 /*--------------------------------------------------------------------------*/
44
45 /* fill_bound_info_3d(): fill boundary info on faces from MACRO_DATA
46 * *data into corresponding macro elements, set boundary types of
47 * vertices (edges are handled later in fill_more_bound_info_3d). The
48 * boundary type of a vertex or edge is just a bit-field where bits
49 * are set according to the boundary types of the surrounding faces.
50 */
51 static void
fill_bound_info_3d(MESH * mesh,const int * mel_vertices,int nv,int ne)52 fill_bound_info_3d(MESH *mesh, const int *mel_vertices, int nv, int ne)
53 {
54 MACRO_EL *mel = mesh->macro_els;
55 int i, j, k;
56 BNDRY_FLAGS bound[nv], non_per_bound[nv];
57 BNDRY_TYPE m_bound;
58
59 /*--------------------------------------------------------------------------*/
60 /* The first N_FACES_3D entries of mel[].boundary will contain pointers to */
61 /* BOUNDARY structures corresponding to the elements' faces. These are the */
62 /* crucial infos provided externally, the other subsimplices (here vertices)*/
63 /* are set to a boundary type depending on the face containing them. */
64 /*--------------------------------------------------------------------------*/
65
66 for (i = 0; i < nv; i++) {
67 BNDRY_FLAGS_INIT(non_per_bound[i]);
68 BNDRY_FLAGS_INIT(bound[i]);
69 }
70
71 for (i = 0; i < ne; i++) {
72 for (j = 0; j < N_WALLS_3D; j++) {
73 if ((m_bound = mel[i].wall_bound[j]) != INTERIOR) {
74 if (mel[i].neigh_vertices[j][0] == -1) {
75 for (k = 1; k < N_VERTICES_3D; k++) {
76 BNDRY_FLAGS_SET(bound[mel_vertices[VERT_IND(3, i, (j+k)%4)]],
77 m_bound);
78 }
79 }
80 for (k = 1; k < N_VERTICES_3D; k++) {
81 BNDRY_FLAGS_SET(non_per_bound[mel_vertices[VERT_IND(3, i, (j+k)%4)]],
82 m_bound);
83 }
84 }
85 }
86 }
87
88 for (i = 0; i < ne; i++) {
89 for (j = 0; j < N_VERTICES_3D; j++) {
90 BNDRY_FLAGS_CPY(mel[i].vertex_bound[j],
91 bound[mel_vertices[VERT_IND(3, i, j)]]);
92 BNDRY_FLAGS_CPY(mel[i].np_vertex_bound[j],
93 non_per_bound[mel_vertices[VERT_IND(3, i, j)]]);
94 }
95 }
96 }
97
98
99 /*--------------------------------------------------------------------------*/
100 /* new_edge_3d(): This procedure counts the number of faces */
101 /* sharing the edge. Another important task is setting the boundary type of */
102 /* the edge depending on the types of these faces via new_bound_val_3d(). */
103 /*--------------------------------------------------------------------------*/
104
105 static inline
new_edge_3d(const int * mel_vertices,MACRO_EL * mel,int mel_edge_no,int * n_neigh,int non_periodic)106 int new_edge_3d(const int *mel_vertices,
107 MACRO_EL *mel, int mel_edge_no, int *n_neigh,
108 int non_periodic)
109 {
110 FUNCNAME("new_edge_3d");
111 typedef BNDRY_FLAGS *LIST_BOUND;
112 MACRO_EL *neigh, *old_neigh = NULL;
113 int i, j, k, opp_v = -1, mel_index, edge_no;
114 int vertex[2] = { -1, }, face, dir;
115 BNDRY_FLAGS bound;
116 int max_no_list_el = 20;
117 LIST_BOUND *list_bound = MEM_ALLOC(max_no_list_el, LIST_BOUND);
118 int np_n_neigh, neigh_periodic = false;
119 static const int next_el[6][2] = {{3,2},{1,3},{1,2},{0,3},{0,2},{0,1}};
120
121 if (non_periodic) {
122 n_neigh = &np_n_neigh;
123 }
124
125 *n_neigh = 1;
126
127 mel_index = mel->index;
128
129 BNDRY_FLAGS_INIT(bound);
130 if (non_periodic) {
131 list_bound[0] = &(mel->np_edge_bound[mel_edge_no]);
132 } else {
133 list_bound[0] = &(mel->edge_bound[mel_edge_no]);
134 }
135
136 /*--------------------------------------------------------------------------*/
137 /* first look for all neighbours in one direction until a boundary is */
138 /* reached :-( or we are back to mel :-) */
139 /* if the index of a neighbour is smaller than the element index, the edge */
140 /* is counted by this neighbour, return 0. */
141 /* If we are back to element, return 1, to count the edge */
142 /*--------------------------------------------------------------------------*/
143
144 neigh = NULL;
145 for (dir = 0; dir < 2; dir++) {
146
147 if (!neigh) {
148
149 for (j = 0; j < 2; j++) {
150 vertex[j] = mel_vertices[VERT_IND(3, mel_index,
151 vertex_of_edge_3d[mel_edge_no][j])];
152 }
153
154 edge_no = mel_edge_no;
155 face = next_el[edge_no][dir];
156 neigh = mel->neigh[face];
157 neigh_periodic = mel->neigh_vertices[face][0] != -1;
158 if (non_periodic && neigh_periodic) {
159 neigh = NULL; /* neighbour across periodic wall */
160 }
161
162 opp_v = mel->opp_vertex[face];
163 if ((non_periodic || !neigh_periodic) && mel->wall_bound[face] != INTERIOR) {
164 BNDRY_FLAGS_SET(bound, mel->wall_bound[face]);
165 }
166
167 old_neigh = mel;
168 }
169
170 while (neigh && neigh != mel) {
171 if (neigh->index < mel_index) {
172 MEM_FREE(list_bound, max_no_list_el, LIST_BOUND);
173 return false;
174 }
175
176 if (neigh_periodic) {
177 for (i = 0; i < N_VERTICES_2D; i++) {
178 j = neigh->neigh_vertices[opp_v][i];
179 if (mel_vertices[VERT_IND(3, old_neigh->index, j)] == vertex[0])
180 break;
181 }
182 j = (opp_v + i + 1) % N_VERTICES_3D;
183 vertex[0] = mel_vertices[VERT_IND(3, neigh->index, j)];
184 for (i = 0; i < N_VERTICES_2D; i++) {
185 k = neigh->neigh_vertices[opp_v][i];
186 if (mel_vertices[VERT_IND(3, old_neigh->index, k)] == vertex[1])
187 break;
188 }
189 k = (opp_v + i + 1) % N_VERTICES_3D;
190 vertex[1] = mel_vertices[VERT_IND(3, neigh->index, k)];
191 } else {
192 for (i = 0; i < N_VERTICES_2D; i++) {
193 j = (opp_v + i + 1) % N_VERTICES_3D;
194 if (mel_vertices[VERT_IND(3, neigh->index, j)] == vertex[0])
195 break;
196 }
197 for (i = 0; i < N_VERTICES_2D; i++) {
198 k = (opp_v + i + 1) % N_VERTICES_3D;
199 if (mel_vertices[VERT_IND(3, neigh->index, k)] == vertex[1])
200 break;
201 }
202 }
203 edge_no = edge_of_vertices_3d[j][k];
204
205 if (*n_neigh == max_no_list_el) {
206 list_bound = MEM_REALLOC(list_bound,
207 max_no_list_el, max_no_list_el + 20,
208 LIST_BOUND);
209 max_no_list_el += 20;
210 }
211
212 if (non_periodic) {
213 list_bound[(*n_neigh)++] = &(neigh->np_edge_bound[edge_no]);
214 } else {
215 list_bound[(*n_neigh)++] = &(neigh->edge_bound[edge_no]);
216 }
217
218 face = (next_el[edge_no][0] == opp_v)
219 ? next_el[edge_no][1] : next_el[edge_no][0];
220
221 neigh_periodic = neigh->neigh_vertices[face][0] != -1;
222 if ((non_periodic || !neigh_periodic) && neigh->wall_bound[face] != INTERIOR) {
223 BNDRY_FLAGS_SET(bound, neigh->wall_bound[face]);
224 }
225
226 opp_v = neigh->opp_vertex[face];
227 if (non_periodic && neigh_periodic) {
228 neigh = NULL; /* neighbour across periodic wall */
229 } else {
230 old_neigh = neigh;
231 neigh = neigh->neigh[face];
232 }
233 }
234 }
235
236 for (j = 0; j < *n_neigh; j++) {
237 BNDRY_FLAGS_CPY(*(list_bound[j]), bound);
238 }
239
240 MEM_FREE(list_bound, max_no_list_el, LIST_BOUND);
241
242 return true;
243 }
244
245 /*--------------------------------------------------------------------------*/
246 /* fill_more_bound_info_3d(): This function fills boundary information using*/
247 /* the function new_edge_3d() above. It also counts faces and edges. */
248 /*--------------------------------------------------------------------------*/
249
250 static void
fill_more_bound_info_3d(MESH * mesh,const int * mel_vertices,int ne,bool do_count)251 fill_more_bound_info_3d(MESH *mesh,
252 const int *mel_vertices, int ne, bool do_count)
253 {
254 /* FUNCNAME("fill_more_bound_info_3d"); */
255 int i, k;
256 int n_edges = 0, n_faces = 0, max_n_neigh = 0, n_neigh = 0;
257 MACRO_EL *mel = mesh->macro_els, *neigh;
258
259 if (do_count) {
260 mesh->max_edge_neigh = 8;
261 }
262
263 for (i = 0; i < ne; i++) {
264 for (k = 0; k < N_EDGES_3D; k++) {
265 /*--------------------------------------------------------------------------*/
266 /* check for not counted edges */
267 /*--------------------------------------------------------------------------*/
268 if (new_edge_3d(mel_vertices, mel+i, k, &n_neigh, false)) {
269 n_edges++;
270 max_n_neigh = MAX(max_n_neigh, n_neigh);
271 }
272 }
273
274 for (k = 0; k < N_NEIGH_3D; k++) {
275 neigh = mel[i].neigh[k];
276 /*--------------------------------------------------------------------------*/
277 /* face is counted by the element with bigger index */
278 /*--------------------------------------------------------------------------*/
279 if (neigh && (neigh->index > mel[i].index))
280 continue;
281
282 n_faces++;
283 }
284 }
285
286 if (mesh->is_periodic) {
287 /* fill also mel->np_edge_bound for non-periodic mesh-traversal,
288 * and update mesh->n_edges to include edges across periodic
289 * walls. The reason is that n_edges is used to layout storage for
290 * DOF-pointers and DOFs are duplicated across periodic walls.
291 */
292 if (do_count) {
293 mesh->per_n_edges = n_edges;
294 mesh->per_n_faces = n_faces;
295 }
296
297 /* per_n_vertices is computed elsewhere */
298
299 n_edges = 0;
300 n_faces = 0;
301 for (i = 0; i < ne; i++) {
302 for (k = 0; k < N_EDGES_3D; k++) {
303 if (new_edge_3d(mel_vertices, mel+i, k, &n_neigh, true)) {
304 n_edges++;
305 max_n_neigh = MAX(max_n_neigh, n_neigh);
306 }
307 }
308 for (k = 0; k < N_NEIGH_3D; k++) {
309 /* skip periodic neighbours this time; this way n_faces always
310 * corresponds to the number of face DOF-pointers.
311 */
312 neigh = mel[i].neigh_vertices[k][0] == -1 ? mel[i].neigh[k] : NULL;
313
314 /* face is counted by the element with bigger index */
315 if (neigh && (neigh->index > mel[i].index))
316 continue;
317
318 n_faces++;
319 }
320 }
321 }
322
323 if (do_count) {
324 mesh->max_edge_neigh = MAX(mesh->max_edge_neigh, 2*max_n_neigh);
325 mesh->n_edges = n_edges;
326 mesh->n_faces = n_faces;
327 }
328
329 return;
330 }
331
332 /*--------------------------------------------------------------------------*/
333 /* AI_get_orientation_3d(): This function calculates an elements orientation*/
334 /* and returns it as a signed integer suitable for the MACRO_EL orientation */
335 /* field. Also called from read_mesh.c */
336 /*--------------------------------------------------------------------------*/
337
AI_get_orientation_3d(MACRO_EL * mel)338 S_CHAR AI_get_orientation_3d(MACRO_EL *mel)
339 {
340 int i;
341 REAL_DD a;
342
343 if (DIM_OF_WORLD != 3) {
344 return 1;
345 }
346
347 for (i = 0; i < 3; i++) {
348 AXPBY_DOW(1.0, *mel->coord[i+1], -1.0, *mel->coord[0], a[i]);
349 }
350
351 return MDET_DOW((const REAL_D *)a) >= 0 ? 1 : -1;
352 }
353
354 /*--------------------------------------------------------------------------*/
355 /* macro_test_3d(): Check data for potential cycles during refinement */
356 /* At the moment, a correction (and subsequent writing of a correct macro */
357 /* data file) can only be done in 2D. */
358 /*--------------------------------------------------------------------------*/
359
360 #if 0
361 ALBERTA_UNUSED(static int check_cycles_3d(MACRO_DATA *));
362 #endif
363
macro_test_3d(MACRO_DATA * data,const char * new_name)364 static void macro_test_3d(MACRO_DATA *data, const char *new_name)
365 {
366 FUNCNAME("macro_test_3d");
367 #if 1
368 WARNING("not implemented for 3d yet: no check is performed\n");
369 #else
370 U_CHAR error_found = false;
371 int i = -1;
372
373 i = check_cycles_3d(data);
374
375 if (i >= 0)
376 {
377 error_found = true;
378 WARNING("There is a cycle beginning in macro element %d.\n", i);
379 ERROR_EXIT("Correction of cycles not implemented yet in 3D\n");
380 }
381
382 if (error_found && new_name) {
383 MSG("Attempting to write corrected macro data to file %s...\n", new_name);
384 write_macro_data(data, new_name);
385 }
386 #endif
387 return;
388 }
389