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: Common/coarsen.c */
7 /* */
8 /* */
9 /* description: coarsening of ? dim. hierarchical meshes; */
10 /* file includes all coarsen_?d.c files !! */
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 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
28 /* */
29 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
30 /* */
31 /*--------------------------------------------------------------------------*/
32
33 #include <stdlib.h>
34 #include "alberta_intern.h"
35 #include "alberta.h"
36
37 /*--------------------------------------------------------------------------*/
38 /* restriction of dof vectors during mesh coarsening */
39 /* using DOF_VEC_LIST structure */
40 /*--------------------------------------------------------------------------*/
41
42 static int
count_coarse_restrict(MESH * mesh,DOF_VEC_LIST * dvlist,int non_periodic)43 count_coarse_restrict(MESH *mesh, DOF_VEC_LIST *dvlist, int non_periodic)
44 {
45 FUNCNAME("count_coarse_restrict");
46 DOF_ADMIN *admin;
47 int iadmin;
48 DOF_REAL_VEC *drv;
49 DOF_REAL_D_VEC *drdv;
50 DOF_INT_VEC *div;
51 DOF_DOF_VEC *ddv;
52 DOF_UCHAR_VEC *duv;
53 DOF_SCHAR_VEC *dsv;
54 DOF_PTR_VEC *dpv;
55 DOF_MATRIX *dm;
56 int ncr = 0;
57 int n_dof_int_vec = 0, n_dof_dof_vec = 0;
58 int n_dof_uchar_vec = 0, n_dof_schar_vec = 0;
59 int n_dof_real_vec = 0, n_dof_real_d_vec = 0;
60 int n_dof_ptr_vec = 0;
61 int n_dof_matrix = 0, n_dof_dowb_matrix = 0, n_dof_rdr_matrix = 0;
62
63 for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
64 admin = mesh->dof_admin[iadmin];
65
66 if (mesh->is_periodic) {
67 if (non_periodic && (admin->flags & ADM_PERIODIC)) {
68 continue;
69 } else if (!non_periodic && !(admin->flags & ADM_PERIODIC)) {
70 continue;
71 }
72 }
73
74 for (div = admin->dof_int_vec; div; div = div->next) {
75 if (div->coarse_restrict) n_dof_int_vec++;
76 }
77 for (ddv = admin->dof_dof_vec; ddv; ddv = ddv->next) {
78 if (ddv->coarse_restrict) n_dof_dof_vec++;
79 }
80 for (ddv = admin->int_dof_vec; ddv; ddv = ddv->next) {
81 if (ddv->coarse_restrict) n_dof_int_vec++;
82 }
83
84 for (duv = admin->dof_uchar_vec; duv; duv = duv->next) {
85 if (duv->coarse_restrict) n_dof_uchar_vec++;
86 }
87 for (dsv = admin->dof_schar_vec; dsv; dsv = dsv->next) {
88 if (dsv->coarse_restrict) n_dof_schar_vec++;
89 }
90
91 for (drv = admin->dof_real_vec; drv; drv = drv->next) {
92 if (drv->coarse_restrict) n_dof_real_vec++;
93 }
94 for (drdv = admin->dof_real_d_vec; drdv; drdv = drdv->next) {
95 if (drdv->coarse_restrict) n_dof_real_d_vec++;
96 }
97
98 for (dpv = admin->dof_ptr_vec; dpv; dpv = dpv->next) {
99 if (dpv->coarse_restrict) n_dof_ptr_vec++;
100 }
101
102 for (dm = admin->dof_matrix; dm; dm = dm->next) {
103 if (dm->coarse_restrict) n_dof_matrix++;
104 }
105 }
106
107 ncr =
108 n_dof_int_vec + n_dof_dof_vec + n_dof_uchar_vec + n_dof_schar_vec
109 + n_dof_ptr_vec + n_dof_real_d_vec
110 + n_dof_real_vec
111 + n_dof_matrix + n_dof_rdr_matrix + n_dof_dowb_matrix;
112
113 if (ncr > 0) {
114 if (dvlist->size < ncr) {
115 dvlist->list = MEM_REALLOC(dvlist->list, dvlist->size, ncr+5, void *);
116 dvlist->size = ncr+5;
117 }
118
119 ncr = 0;
120
121 if (n_dof_int_vec)
122 dvlist->dof_int_vec = (DOF_INT_VEC **)(dvlist->list+ncr);
123 else
124 dvlist->dof_int_vec = NULL;
125 ncr += n_dof_int_vec;
126
127 if (n_dof_dof_vec)
128 dvlist->dof_dof_vec = (DOF_DOF_VEC **)(dvlist->list+ncr);
129 else
130 dvlist->dof_dof_vec = NULL;
131 ncr += n_dof_dof_vec;
132
133 if (n_dof_uchar_vec)
134 dvlist->dof_uchar_vec = (DOF_UCHAR_VEC **)(dvlist->list+ncr);
135 else
136 dvlist->dof_uchar_vec = NULL;
137 ncr += n_dof_uchar_vec;
138
139 if (n_dof_schar_vec)
140 dvlist->dof_schar_vec = (DOF_SCHAR_VEC **)(dvlist->list+ncr);
141 else
142 dvlist->dof_schar_vec = NULL;
143 ncr += n_dof_schar_vec;
144
145 if (n_dof_real_vec)
146 dvlist->dof_real_vec = (DOF_REAL_VEC **)(dvlist->list+ncr);
147 else
148 dvlist->dof_real_vec = NULL;
149 ncr += n_dof_real_vec;
150
151 if (n_dof_real_d_vec)
152 dvlist->dof_real_d_vec = (DOF_REAL_D_VEC **)(dvlist->list+ncr);
153 else
154 dvlist->dof_real_d_vec = NULL;
155 ncr += n_dof_real_d_vec;
156
157 if (n_dof_ptr_vec)
158 dvlist->dof_ptr_vec = (DOF_PTR_VEC **)(dvlist->list+ncr);
159 else
160 dvlist->dof_ptr_vec = NULL;
161 ncr += n_dof_ptr_vec;
162
163 if (n_dof_matrix)
164 dvlist->dof_matrix = (DOF_MATRIX **)(dvlist->list+ncr);
165 else
166 dvlist->dof_matrix = NULL;
167 ncr += n_dof_matrix;
168
169 DEBUG_TEST_EXIT(ncr <= dvlist->size, "error in dvlist->size");
170
171 dvlist->n_dof_int_vec = 0;
172 dvlist->n_dof_dof_vec = 0;
173 dvlist->n_dof_uchar_vec = 0;
174 dvlist->n_dof_schar_vec = 0;
175 dvlist->n_dof_real_vec = 0;
176 dvlist->n_dof_real_d_vec = 0;
177 dvlist->n_dof_ptr_vec = 0;
178 dvlist->n_dof_matrix = 0;
179
180 for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
181 admin = mesh->dof_admin[iadmin];
182
183 if (mesh->is_periodic) {
184 if (non_periodic && (admin->flags & ADM_PERIODIC)) {
185 continue;
186 } else if (!non_periodic && !(admin->flags & ADM_PERIODIC)) {
187 continue;
188 }
189 }
190
191 for (div = admin->dof_int_vec; div; div = div->next) {
192 if (div->coarse_restrict)
193 dvlist->dof_int_vec[dvlist->n_dof_int_vec++] = div;
194 }
195 for (ddv = admin->dof_dof_vec; ddv; ddv = ddv->next) {
196 if (ddv->coarse_restrict)
197 dvlist->dof_dof_vec[dvlist->n_dof_dof_vec++] = ddv;
198 }
199 for (ddv = admin->int_dof_vec; ddv; ddv = ddv->next) {
200 if (ddv->coarse_restrict)
201 dvlist->dof_dof_vec[dvlist->n_dof_dof_vec++] = ddv;
202 }
203 for (duv = admin->dof_uchar_vec; duv; duv = duv->next) {
204 if (duv->coarse_restrict)
205 dvlist->dof_uchar_vec[dvlist->n_dof_uchar_vec++] = duv;
206 }
207 for (dsv = admin->dof_schar_vec; dsv; dsv = dsv->next) {
208 if (dsv->coarse_restrict)
209 dvlist->dof_schar_vec[dvlist->n_dof_schar_vec++] = dsv;
210 }
211 for (drv = admin->dof_real_vec; drv; drv = drv->next) {
212 if (drv->coarse_restrict)
213 dvlist->dof_real_vec[dvlist->n_dof_real_vec++] = drv;
214 }
215 for (drdv = admin->dof_real_d_vec; drdv; drdv = drdv->next) {
216 if (drdv->coarse_restrict)
217 dvlist->dof_real_d_vec[dvlist->n_dof_real_d_vec++] = drdv;
218 }
219 for (dpv = admin->dof_ptr_vec; dpv; dpv = dpv->next) {
220 if (dpv->coarse_restrict)
221 dvlist->dof_ptr_vec[dvlist->n_dof_ptr_vec++] = dpv;
222 }
223 for (dm = admin->dof_matrix; dm; dm = dm->next) {
224 if (dm->coarse_restrict)
225 dvlist->dof_matrix[dvlist->n_dof_matrix++] = dm;
226 }
227 }
228
229 DEBUG_TEST_EXIT(dvlist->n_dof_int_vec == n_dof_int_vec,
230 "error in n_dof_int_vec");
231 DEBUG_TEST_EXIT(dvlist->n_dof_dof_vec == n_dof_dof_vec,
232 "error in n_dof_dof_vec");
233 DEBUG_TEST_EXIT(dvlist->n_dof_uchar_vec == n_dof_uchar_vec,
234 "error in n_dof_uchar_vec");
235 DEBUG_TEST_EXIT(dvlist->n_dof_schar_vec == n_dof_schar_vec,
236 "error in n_dof_schar_vec");
237 DEBUG_TEST_EXIT(dvlist->n_dof_real_vec == n_dof_real_vec,
238 "error in n_dof_real_vec");
239 DEBUG_TEST_EXIT(dvlist->n_dof_real_d_vec == n_dof_real_d_vec,
240 "error in n_dof_real_d_vec");
241 DEBUG_TEST_EXIT(dvlist->n_dof_ptr_vec == n_dof_ptr_vec,
242 "error in n_dof_ptr_vec");
243 DEBUG_TEST_EXIT(dvlist->n_dof_matrix == n_dof_matrix,
244 "error in n_dof_matrix");
245 } else {
246 dvlist->dof_int_vec = NULL;
247 dvlist->dof_dof_vec = NULL;
248 dvlist->dof_uchar_vec = NULL;
249 dvlist->dof_schar_vec = NULL;
250 dvlist->dof_real_vec = NULL;
251 dvlist->dof_real_d_vec = NULL;
252 dvlist->dof_ptr_vec = NULL;
253 dvlist->dof_matrix = NULL;
254 }
255
256 return ncr;
257 }
258
coarse_restrict(MESH * mesh,DOF_VEC_LIST * dvlist,RC_LIST_EL * list,int n_el)259 static void coarse_restrict(MESH *mesh, DOF_VEC_LIST *dvlist,
260 RC_LIST_EL *list, int n_el)
261 {
262 DOF_REAL_VEC *drv;
263 DOF_REAL_D_VEC *drdv;
264 DOF_INT_VEC *div;
265 DOF_DOF_VEC *ddv;
266 DOF_UCHAR_VEC *duv;
267 DOF_SCHAR_VEC *dsv;
268 DOF_PTR_VEC *dpv;
269 DOF_MATRIX *dm;
270 int i;
271
272 DEBUG_TEST_EXIT(((MESH_MEM_INFO *)mesh->mem_info)->dvlist,
273 "Huh? No dvlist found?\n");
274
275 for (i = 0; i < dvlist->n_dof_ptr_vec; i++) {
276 dpv = dvlist->dof_ptr_vec[i];
277 if (dpv->coarse_restrict) dpv->coarse_restrict(dpv, list, n_el);
278 }
279
280 for (i = 0; i < dvlist->n_dof_int_vec; i++) {
281 div = dvlist->dof_int_vec[i];
282 if (div->coarse_restrict) div->coarse_restrict(div, list, n_el);
283 }
284 for (i = 0; i < dvlist->n_dof_dof_vec; i++) {
285 ddv = dvlist->dof_dof_vec[i];
286 if (ddv->coarse_restrict) ddv->coarse_restrict(ddv, list, n_el);
287 }
288
289 for (i = 0; i < dvlist->n_dof_uchar_vec; i++) {
290 duv = dvlist->dof_uchar_vec[i];
291 if (duv->coarse_restrict) duv->coarse_restrict(duv, list, n_el);
292 }
293 for (i = 0; i < dvlist->n_dof_schar_vec; i++) {
294 dsv = dvlist->dof_schar_vec[i];
295 if (dsv->coarse_restrict) dsv->coarse_restrict(dsv, list, n_el);
296 }
297
298 for (i = 0; i < dvlist->n_dof_real_vec; i++) {
299 drv = dvlist->dof_real_vec[i];
300 if (drv->coarse_restrict) drv->coarse_restrict(drv, list, n_el);
301 }
302 for (i = 0; i < dvlist->n_dof_real_d_vec; i++) {
303 drdv = dvlist->dof_real_d_vec[i];
304 if (drdv->coarse_restrict) drdv->coarse_restrict(drdv, list, n_el);
305 }
306
307 for (i = 0; i < dvlist->n_dof_matrix; i++) {
308 dm = dvlist->dof_matrix[i];
309 if (dm->coarse_restrict) dm->coarse_restrict(dm, list, n_el);
310 }
311 }
312
313 /*--------------------------------------------------------------------------*/
314 /* calculate the maximal level of the triangulation */
315 /*--------------------------------------------------------------------------*/
316
317 typedef struct coarsen_traverse_data {
318 int max_level_val;
319 U_CHAR global_mark;
320 } COARSEN_TRAVERSE_DATA;
321
322
max_level_fct(const EL_INFO * el_info,void * data)323 static void max_level_fct(const EL_INFO *el_info, void *data)
324 {
325 ((COARSEN_TRAVERSE_DATA *)data)->max_level_val =
326 MAX(((COARSEN_TRAVERSE_DATA *)data)->max_level_val, el_info->level);
327 }
328
get_max_level(MESH * mesh)329 int get_max_level(MESH *mesh)
330 {
331 COARSEN_TRAVERSE_DATA td[1] = {{0}};
332
333 mesh_traverse(mesh, -1, CALL_LEAF_EL, max_level_fct, td);
334 return td->max_level_val;
335 }
336
337 /*--------------------------------------------------------------------------*/
338 /* global_coarsen: */
339 /* tries to coarsen every element of mesh at least mark times */
340 /*--------------------------------------------------------------------------*/
341
global_coarsen(MESH * mesh,int mark,FLAGS fill_flags)342 U_CHAR global_coarsen(MESH *mesh, int mark, FLAGS fill_flags)
343 {
344 if (mark >= 0) {
345 return 0;
346 }
347 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) {
348 el_info->el->mark = mark;
349 } TRAVERSE_NEXT();
350 return coarsen(mesh, fill_flags);
351 }
352
353 static int call_coarse_restrict_1d;
354 #include "coarsen_1d.c"
355
356 #if DIM_MAX > 1
357 static int call_coarse_restrict_2d;
358 static int call_coarse_restrict_np_2d;
359
360 static int do_more_coarsen_2d = true;
361
362 #include "coarsen_2d.c"
363
364 #if DIM_MAX > 2
365 static int call_coarse_restrict_3d;
366 static int call_coarse_restrict_np_3d;
367
368 static int do_more_coarsen_3d = true;
369 static TRAVERSE_STACK *stack_3d;
370 #include "coarsen_3d.c"
371 #endif
372 #endif
373
374
transfer_fct(const EL_INFO * elinfo,void * data)375 static void transfer_fct(const EL_INFO *elinfo, void *data)
376 {
377 MESH *mesh = elinfo->mesh;
378 EL *s_el = elinfo->el,
379 *m_el;
380 DOF_PTR_VEC *master_binding;
381 const DOF_ADMIN *s_admin;
382
383 if(s_el->mark < 0) {
384 master_binding = ((MESH_MEM_INFO *)mesh->mem_info)->master_binding;
385 s_admin = master_binding->fe_space->admin;
386
387 m_el = (EL *)master_binding->vec[s_el->dof[mesh->node[CENTER]]
388 [s_admin->n0_dof[CENTER]]];
389
390 /****************************************************************************/
391 /* Transfer the coarsening marks to the master mesh. PLEASE NOTE: This */
392 /* actually overwrites master mesh marks! Otherwise, we could never do */
393 /* any coarsening for slave meshes. This may not be what the user wishes */
394 /* and should be well documented! */
395 /****************************************************************************/
396 m_el->mark = MIN(m_el->mark, -1);
397 }
398 return;
399 }
400
blank_marks_fct(const EL_INFO * elinfo,void * data)401 static void blank_marks_fct(const EL_INFO *elinfo, void *data)
402 {
403 EL *s_el = elinfo->el;
404
405 s_el->mark = MAX(s_el->mark, 0);
406 return;
407 }
408
coarsen(MESH * mesh,FLAGS fill_flags)409 U_CHAR coarsen(MESH *mesh, FLAGS fill_flags)
410 {
411 FUNCNAME("coarsen");
412 MESH_MEM_INFO *mem_info = ((MESH_MEM_INFO *)mesh->mem_info);
413 U_CHAR mesh_coarsened = 0;
414
415 if (mem_info->n_slaves) {
416 /****************************************************************************/
417 /* We are on a master mesh. */
418 /****************************************************************************/
419
420 #if DIM_MAX > 1
421 /****************************************************************************/
422 /* Check if we have an entire hierarchy of meshes. In this case, we need */
423 /* to set call_coarse_restrict_2d to nonzero to make the triple coarsening */
424 /* work! */
425 /****************************************************************************/
426 if (mesh->dim == 2) {
427 int i;
428 call_coarse_restrict_1d = false;
429
430 for(i = 0; i < mem_info->n_slaves; i++) {
431 MESH *slave = mem_info->slaves[i];
432 call_coarse_restrict_1d +=
433 count_coarse_restrict(slave, AI_get_dof_vec_list(slave), false);
434 }
435 }
436
437 #if DIM_MAX == 3
438 if (mesh->dim == 3) {
439 int i, j;
440 call_coarse_restrict_1d = 0;
441 call_coarse_restrict_2d = 0;
442
443 for(i = 0; i < mem_info->n_slaves; i++) {
444 MESH *slave = mem_info->slaves[i];
445 MESH_MEM_INFO *mem_info_2 = (MESH_MEM_INFO *)slave->mem_info;
446
447 call_coarse_restrict_2d +=
448 count_coarse_restrict(slave, AI_get_dof_vec_list(slave), false);
449 if (slave->is_periodic)
450 call_coarse_restrict_np_2d +=
451 count_coarse_restrict(slave, AI_get_dof_vec_list_np(slave), true);
452
453 for(j = 0; j < mem_info_2->n_slaves; j++) {
454 MESH *slave_2 = mem_info_2->slaves[j];
455
456 call_coarse_restrict_1d +=
457 count_coarse_restrict(slave_2, AI_get_dof_vec_list(slave_2), false);
458 }
459 }
460 }
461 #endif
462 #endif
463 }
464
465 if (mem_info->master) {
466 /****************************************************************************/
467 /* We are on a slave mesh. */
468 /****************************************************************************/
469 int n_slave_elements = mesh->n_elements;
470
471 /* Transfer the coarsening marks to the master mesh. See note above! */
472 do {
473 mesh_traverse(mesh, 0, FILL_NOTHING|CALL_LEAF_EL, transfer_fct, NULL);
474 mesh_coarsened = coarsen(mem_info->master, fill_flags);
475 } while (mesh_coarsened);
476
477 /* Delete the unfortunate remaining coarsening marks. */
478 mesh_traverse(mesh, 0, FILL_NOTHING|CALL_LEAF_EL, blank_marks_fct, NULL);
479
480 return (mesh->n_elements < n_slave_elements) ? MESH_COARSENED : 0;
481 #if 0
482 } else {
483 /****************************************************************************/
484 /* We are on a top level master mesh. */
485 /****************************************************************************/
486
487 /* Advance cookies on this mesh and all its slaves. */
488 AI_advance_cookies_rec(mesh);
489 #endif
490 }
491
492 switch (mesh->dim) {
493 case 0:
494 WARNING("No coarsening possible for dim==0!\n");
495 break;
496 case 1:
497 return coarsen_1d(mesh, fill_flags);
498 break;
499 #if DIM_MAX > 1
500 case 2:
501 return coarsen_2d(mesh, fill_flags);
502 break;
503 #if DIM_MAX > 2
504 case 3:
505 return coarsen_3d(mesh, fill_flags);
506 break;
507 #endif
508 #endif
509 default:
510 ERROR_EXIT("Illegal dim during coarsening!\n");
511 break;
512 }
513
514 if (mesh_coarsened) {
515 /* Advance cookies on this mesh and all its slaves. */
516 AI_advance_cookies_rec(mesh);
517 }
518
519 return mesh_coarsened;
520 }
521