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/refine.c */
7 /* */
8 /* */
9 /* description: Common refinement routines shared among all dimensions */
10 /* This file contains refine_?d.c for ? = 1,2,3 !! */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Alfred Schmidt */
14 /* Zentrum fuer Technomathematik */
15 /* Fachbereich 3 Mathematik/Informatik */
16 /* Univesitaet 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 #include "alberta_intern.h"
33 #include "alberta.h"
34
35 /*--------------------------------------------------------------------------*/
36 /* refinement routines for global refinement: */
37 /*--------------------------------------------------------------------------*/
38
39 /*--------------------------------------------------------------------------*/
40 /* global_refine: */
41 /* refines every element of mesh at least mark times */
42 /*--------------------------------------------------------------------------*/
43
global_refine(MESH * mesh,int mark,FLAGS fill_flags)44 U_CHAR global_refine(MESH *mesh, int mark, FLAGS fill_flags)
45 {
46 if (mark <= 0) {
47 return 0;
48 }
49 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) {
50 el_info->el->mark = mark;
51 } TRAVERSE_NEXT();
52
53 return refine(mesh, fill_flags);
54 }
55
56
AI_get_dof_vec_list(MESH * mesh)57 DOF_VEC_LIST *AI_get_dof_vec_list(MESH *mesh)
58 {
59 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
60
61 if (!mem_info->dvlist)
62 mem_info->dvlist = MEM_CALLOC(1, DOF_VEC_LIST);
63
64 return mem_info->dvlist;
65 }
66
AI_free_dof_vec_list(MESH * mesh)67 void AI_free_dof_vec_list(MESH *mesh)
68 {
69 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
70
71 if (mem_info->dvlist) {
72 MEM_FREE(mem_info->dvlist, 1, DOF_VEC_LIST);
73 mem_info->dvlist = NULL;
74 }
75 }
76
AI_get_dof_vec_list_np(MESH * mesh)77 DOF_VEC_LIST *AI_get_dof_vec_list_np(MESH *mesh)
78 {
79 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
80
81 if (!mem_info->dvlist_np)
82 mem_info->dvlist_np = MEM_CALLOC(1, DOF_VEC_LIST);
83
84 return mem_info->dvlist_np;
85 }
86
AI_free_dof_vec_list_np(MESH * mesh)87 void AI_free_dof_vec_list_np(MESH *mesh)
88 {
89 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
90
91 if (mem_info->dvlist_np) {
92 MEM_FREE(mem_info->dvlist_np, 1, DOF_VEC_LIST);
93 mem_info->dvlist_np = NULL;
94 }
95 }
96
97 /* Determine how many vectors/matrices request an adaptation
98 * call-back. Fill FILL_FLAGS with the fill-flag requirements of the
99 * underlying basis functions. FILL_FLAGS is assumed to be initialized
100 * by the caller.
101 */
102 static int
count_refine_interpol(MESH * mesh,DOF_VEC_LIST * dvlist,int non_periodic,FLAGS * fill_flags)103 count_refine_interpol(MESH *mesh, DOF_VEC_LIST *dvlist,
104 int non_periodic, FLAGS *fill_flags)
105 {
106 FUNCNAME("count_refine_interpol");
107 DOF_ADMIN *admin;
108 int iadmin;
109 DOF_REAL_VEC *drv;
110 DOF_REAL_D_VEC *drdv;
111 DOF_INT_VEC *div;
112 DOF_DOF_VEC *ddv;
113 DOF_UCHAR_VEC *duv;
114 DOF_SCHAR_VEC *dsv;
115 DOF_PTR_VEC *dpv;
116 DOF_MATRIX *dm;
117 int nri = 0;
118 int n_dof_int_vec = 0, n_dof_dof_vec = 0;
119 int n_dof_uchar_vec = 0, n_dof_schar_vec = 0;
120 int n_dof_real_vec = 0, n_dof_real_d_vec = 0;
121 int n_dof_ptr_vec = 0;
122 int n_dof_matrix = 0, n_dof_dowb_matrix = 0, n_dof_rdr_matrix = 0;
123
124 for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
125 admin = mesh->dof_admin[iadmin];
126
127 if (mesh->is_periodic) {
128 if (non_periodic && (admin->flags & ADM_PERIODIC)) {
129 continue;
130 } else if (!non_periodic && !(admin->flags & ADM_PERIODIC)) {
131 continue;
132 }
133 }
134
135 for (div = admin->dof_int_vec; div; div = div->next) {
136 if (div->refine_interpol) {
137 n_dof_int_vec++;
138 if (div->fe_space->bas_fcts) {
139 *fill_flags |= div->fe_space->bas_fcts->fill_flags;
140 }
141 }
142 }
143 for (ddv = admin->dof_dof_vec; ddv; ddv = ddv->next) {
144 if (ddv->refine_interpol) {
145 n_dof_dof_vec++;
146 if (ddv->fe_space->bas_fcts) {
147 *fill_flags |= ddv->fe_space->bas_fcts->fill_flags;
148 }
149 }
150 }
151 for (ddv = admin->int_dof_vec; ddv; ddv = ddv->next) {
152 if (ddv->refine_interpol) {
153 n_dof_dof_vec++;
154 if (ddv->fe_space->bas_fcts) {
155 *fill_flags |= ddv->fe_space->bas_fcts->fill_flags;
156 }
157 }
158 }
159
160 for (duv = admin->dof_uchar_vec; duv; duv = duv->next) {
161 if (duv->refine_interpol) {
162 n_dof_uchar_vec++;
163 if (duv->fe_space->bas_fcts) {
164 *fill_flags |= duv->fe_space->bas_fcts->fill_flags;
165 }
166 }
167 }
168 for (dsv = admin->dof_schar_vec; dsv; dsv = dsv->next) {
169 if (dsv->refine_interpol) {
170 n_dof_schar_vec++;
171 if (dsv->fe_space->bas_fcts) {
172 *fill_flags |= dsv->fe_space->bas_fcts->fill_flags;
173 }
174 }
175 }
176
177 for (drv = admin->dof_real_vec; drv; drv = drv->next) {
178 if (drv->refine_interpol) {
179 n_dof_real_vec++;
180 if (drv->fe_space->bas_fcts) {
181 *fill_flags |= drv->fe_space->bas_fcts->fill_flags;
182 }
183 }
184 }
185 for (drdv = admin->dof_real_d_vec; drdv; drdv = drdv->next) {
186 if (drdv->refine_interpol) {
187 n_dof_real_d_vec++;
188 if (drdv->fe_space->bas_fcts) {
189 *fill_flags |= drdv->fe_space->bas_fcts->fill_flags;
190 }
191 }
192 }
193
194 for (dpv = admin->dof_ptr_vec; dpv; dpv = dpv->next) {
195 if (dpv->refine_interpol) {
196 n_dof_ptr_vec++;
197 if (dpv->fe_space->bas_fcts) {
198 *fill_flags |= dpv->fe_space->bas_fcts->fill_flags;
199 }
200 }
201 }
202
203 for (dm = admin->dof_matrix; dm; dm = dm->next) {
204 if (dm->refine_interpol) {
205 n_dof_matrix++;
206 if (dm->row_fe_space->bas_fcts) {
207 *fill_flags |= dm->row_fe_space->bas_fcts->fill_flags;
208 }
209 if (dm->col_fe_space->bas_fcts) {
210 *fill_flags |= dm->col_fe_space->bas_fcts->fill_flags;
211 }
212 }
213 }
214 }
215
216 nri =
217 + n_dof_int_vec + n_dof_dof_vec + n_dof_uchar_vec + n_dof_schar_vec
218 + n_dof_real_vec + n_dof_real_d_vec
219 + n_dof_ptr_vec
220 + n_dof_matrix + n_dof_rdr_matrix + n_dof_dowb_matrix;
221
222 if (nri > 0) {
223 if (dvlist->size < nri) {
224 dvlist->list = MEM_REALLOC(dvlist->list, dvlist->size, nri+5, void *);
225 dvlist->size = nri+5;
226 }
227
228 nri = 0;
229
230 if (n_dof_int_vec)
231 dvlist->dof_int_vec = (DOF_INT_VEC **)(dvlist->list+nri);
232 else
233 dvlist->dof_int_vec = NULL;
234 nri += n_dof_int_vec;
235
236 if (n_dof_dof_vec)
237 dvlist->dof_dof_vec = (DOF_DOF_VEC **)(dvlist->list+nri);
238 else
239 dvlist->dof_dof_vec = NULL;
240 nri += n_dof_dof_vec;
241
242 if (n_dof_uchar_vec)
243 dvlist->dof_uchar_vec = (DOF_UCHAR_VEC **)(dvlist->list+nri);
244 else
245 dvlist->dof_uchar_vec = NULL;
246 nri += n_dof_uchar_vec;
247
248 if (n_dof_schar_vec)
249 dvlist->dof_schar_vec = (DOF_SCHAR_VEC **)(dvlist->list+nri);
250 else
251 dvlist->dof_schar_vec = NULL;
252 nri += n_dof_schar_vec;
253
254 if (n_dof_real_vec)
255 dvlist->dof_real_vec = (DOF_REAL_VEC **)(dvlist->list+nri);
256 else
257 dvlist->dof_real_vec = NULL;
258 nri += n_dof_real_vec;
259
260 if (n_dof_real_d_vec)
261 dvlist->dof_real_d_vec = (DOF_REAL_D_VEC **)(dvlist->list+nri);
262 else
263 dvlist->dof_real_d_vec = NULL;
264 nri += n_dof_real_d_vec;
265
266 if (n_dof_ptr_vec)
267 dvlist->dof_ptr_vec = (DOF_PTR_VEC **)(dvlist->list+nri);
268 else
269 dvlist->dof_ptr_vec = NULL;
270 nri += n_dof_ptr_vec;
271
272 if (n_dof_matrix)
273 dvlist->dof_matrix = (DOF_MATRIX **)(dvlist->list+nri);
274 else
275 dvlist->dof_matrix = NULL;
276 nri += n_dof_matrix;
277
278 DEBUG_TEST_EXIT(nri <= dvlist->size, "error in dvlist->size");
279
280 dvlist->n_dof_int_vec = 0;
281 dvlist->n_dof_dof_vec = 0;
282 dvlist->n_dof_uchar_vec = 0;
283 dvlist->n_dof_schar_vec = 0;
284 dvlist->n_dof_real_vec = 0;
285 dvlist->n_dof_real_d_vec = 0;
286 dvlist->n_dof_ptr_vec = 0;
287 dvlist->n_dof_matrix = 0;
288
289 for (iadmin = 0; iadmin < mesh->n_dof_admin; iadmin++) {
290 admin = mesh->dof_admin[iadmin];
291
292 if (mesh->is_periodic) {
293 if (non_periodic && (admin->flags & ADM_PERIODIC)) {
294 continue;
295 } else if (!non_periodic && !(admin->flags & ADM_PERIODIC)) {
296 continue;
297 }
298 }
299
300 for (div = admin->dof_int_vec; div; div = div->next) {
301 if (div->refine_interpol)
302 dvlist->dof_int_vec[dvlist->n_dof_int_vec++] = div;
303 }
304 for (ddv = admin->dof_dof_vec; ddv; ddv = ddv->next) {
305 if (ddv->refine_interpol)
306 dvlist->dof_dof_vec[dvlist->n_dof_dof_vec++] = ddv;
307 }
308 for (ddv = admin->int_dof_vec; ddv; ddv = ddv->next) {
309 if (ddv->refine_interpol)
310 dvlist->dof_dof_vec[dvlist->n_dof_dof_vec++] = ddv;
311 }
312
313 for (duv = admin->dof_uchar_vec; duv; duv = duv->next) {
314 if (duv->refine_interpol)
315 dvlist->dof_uchar_vec[dvlist->n_dof_uchar_vec++] = duv;
316 }
317 for (dsv = admin->dof_schar_vec; dsv; dsv = dsv->next) {
318 if (dsv->refine_interpol)
319 dvlist->dof_schar_vec[dvlist->n_dof_schar_vec++] = dsv;
320 }
321
322 for (drv = admin->dof_real_vec; drv; drv = drv->next) {
323 if (drv->refine_interpol)
324 dvlist->dof_real_vec[dvlist->n_dof_real_vec++] = drv;
325 }
326 for (drdv = admin->dof_real_d_vec; drdv; drdv = drdv->next) {
327 if (drdv->refine_interpol)
328 dvlist->dof_real_d_vec[dvlist->n_dof_real_d_vec++] = drdv;
329 }
330
331 for (dpv = admin->dof_ptr_vec; dpv; dpv = dpv->next) {
332 if (dpv->refine_interpol)
333 dvlist->dof_ptr_vec[dvlist->n_dof_ptr_vec++] = dpv;
334 }
335
336 for (dm = admin->dof_matrix; dm; dm = dm->next) {
337 if (dm->refine_interpol)
338 dvlist->dof_matrix[dvlist->n_dof_matrix++] = dm;
339 }
340 }
341
342 DEBUG_TEST_EXIT(dvlist->n_dof_int_vec == n_dof_int_vec,
343 "error in n_dof_int_vec");
344 DEBUG_TEST_EXIT(dvlist->n_dof_dof_vec == n_dof_dof_vec,
345 "error in n_dof_dof_vec");
346 DEBUG_TEST_EXIT(dvlist->n_dof_uchar_vec == n_dof_uchar_vec,
347 "error in n_dof_uchar_vec");
348 DEBUG_TEST_EXIT(dvlist->n_dof_schar_vec == n_dof_schar_vec,
349 "error in n_dof_schar_vec");
350 DEBUG_TEST_EXIT(dvlist->n_dof_real_vec == n_dof_real_vec,
351 "error in n_dof_real_vec");
352 DEBUG_TEST_EXIT(dvlist->n_dof_real_d_vec == n_dof_real_d_vec,
353 "error in n_dof_real_d_vec");
354 DEBUG_TEST_EXIT(dvlist->n_dof_ptr_vec == n_dof_ptr_vec,
355 "error in n_dof_ptr_vec");
356 DEBUG_TEST_EXIT(dvlist->n_dof_matrix == n_dof_matrix,
357 "error in n_dof_matrix");
358 } else {
359 dvlist->dof_int_vec = NULL;
360 dvlist->dof_dof_vec = NULL;
361 dvlist->dof_uchar_vec = NULL;
362 dvlist->dof_schar_vec = NULL;
363 dvlist->dof_real_vec = NULL;
364 dvlist->dof_real_d_vec = NULL;
365 dvlist->dof_ptr_vec = NULL;
366 dvlist->dof_matrix = NULL;
367 }
368
369 return nri;
370 }
371
refine_interpol(MESH * mesh,DOF_VEC_LIST * dvlist,RC_LIST_EL * list,int n_el)372 static void refine_interpol(MESH *mesh, DOF_VEC_LIST *dvlist,
373 RC_LIST_EL *list, int n_el)
374 {
375 DOF_REAL_VEC *drv;
376 DOF_REAL_D_VEC *drdv;
377 DOF_INT_VEC *div;
378 DOF_DOF_VEC *ddv;
379 DOF_UCHAR_VEC *duv;
380 DOF_SCHAR_VEC *dsv;
381 DOF_PTR_VEC *dpv;
382 DOF_MATRIX *dm;
383 int i;
384
385 for (i = 0; i < dvlist->n_dof_int_vec; i++) {
386 div = dvlist->dof_int_vec[i];
387 div->refine_interpol(div, list, n_el);
388 }
389 for (i = 0; i < dvlist->n_dof_dof_vec; i++) {
390 ddv = dvlist->dof_dof_vec[i];
391 ddv->refine_interpol(ddv, list, n_el);
392 }
393
394 for (i = 0; i < dvlist->n_dof_uchar_vec; i++) {
395 duv = dvlist->dof_uchar_vec[i];
396 duv->refine_interpol(duv, list, n_el);
397 }
398 for (i = 0; i < dvlist->n_dof_schar_vec; i++) {
399 dsv = dvlist->dof_schar_vec[i];
400 dsv->refine_interpol(dsv, list, n_el);
401 }
402
403 for (i = 0; i < dvlist->n_dof_real_vec; i++) {
404 drv = dvlist->dof_real_vec[i];
405 drv->refine_interpol(drv, list, n_el);
406 }
407 for (i = 0; i < dvlist->n_dof_real_d_vec; i++) {
408 drdv = dvlist->dof_real_d_vec[i];
409 drdv->refine_interpol(drdv, list, n_el);
410 }
411
412 /* Refine DOF_PTR_VECs _AFTER_ DOF_REAL_D_VECs; otherwise
413 * slave_refine_interpolY() will not work. It is called from
414 * master_interpol() and needs the refined coordinate information
415 * from the master mesh.
416 */
417 for (i = 0; i < dvlist->n_dof_ptr_vec; i++) {
418 dpv = dvlist->dof_ptr_vec[i];
419 dpv->refine_interpol(dpv, list, n_el);
420 }
421
422 for (i = 0; i < dvlist->n_dof_matrix; i++) {
423 dm = dvlist->dof_matrix[i];
424 dm->refine_interpol(dm, list, n_el);
425 }
426
427 return;
428 }
429
_AI_refine_update_bbox(MESH * mesh,const REAL_D new_coord)430 void _AI_refine_update_bbox(MESH *mesh, const REAL_D new_coord)
431 {
432 int i;
433
434 for (i = 0; i < DIM_OF_WORLD; i++) {
435 if (new_coord[i] < mesh->bbox[0][i]) {
436 mesh->bbox[0][i] = new_coord[i];
437 mesh->diam[i] = mesh->bbox[1][i] - mesh->bbox[0][i];
438 } else if (new_coord[i] > mesh->bbox[1][i]) {
439 mesh->bbox[1][i] = new_coord[i];
440 mesh->diam[i] = mesh->bbox[1][i] - mesh->bbox[0][i];
441 }
442 }
443 }
444
445 static int do_more_refine_1d = 0;
446 static int call_refine_interpol_1d = 0;
447
448 #include "refine_1d.c"
449
450 #if DIM_MAX > 1
451 static int do_more_refine_2d = false;
452 static int call_refine_interpol_2d = false;
453 static int call_refine_interpol_np_2d = false;
454
455 #include "refine_2d.c"
456
457 #if DIM_MAX > 2
458 static int do_more_refine_3d = false;
459 static int call_refine_interpol_3d = false;
460 static int call_refine_interpol_np_3d = false;
461
462 #include "refine_3d.c"
463 #endif
464 #endif
465
466 /*--------------------------------------------------------------------------*/
467 /* refine: */
468 /* traversal routine for recursive refinement of a triangulation; basic */
469 /* refinement module; used by all other refinement routines like */
470 /* global_refine() */
471 /*--------------------------------------------------------------------------*/
472
transfer_fct(const EL_INFO * elinfo,void * data)473 static void transfer_fct(const EL_INFO *elinfo, void *data)
474 {
475 MESH *mesh = elinfo->mesh;
476 EL *s_el = elinfo->el,
477 *m_el;
478 DOF_PTR_VEC *master_binding;
479 const DOF_ADMIN *s_admin;
480
481 if (s_el->mark > 0) {
482 master_binding = ((MESH_MEM_INFO *)mesh->mem_info)->master_binding;
483 s_admin = master_binding->fe_space->admin;
484
485 m_el = (EL *)master_binding->vec[s_el->dof[mesh->node[CENTER]]
486 [s_admin->n0_dof[CENTER]]];
487
488 m_el->mark = MAX(m_el->mark, 1);
489 }
490
491 return;
492 }
493
refine(MESH * mesh,FLAGS fill_flags)494 U_CHAR refine(MESH *mesh, FLAGS fill_flags)
495 {
496 FUNCNAME("refine");
497 MESH_MEM_INFO *mem_info = ((MESH_MEM_INFO *)mesh->mem_info);
498 U_CHAR mesh_refined = 0;
499
500 if (mem_info->n_slaves) {
501 /****************************************************************************/
502 /* We are on a master mesh. */
503 /****************************************************************************/
504
505 #if DIM_MAX > 1
506 /****************************************************************************/
507 /* Check if we have an entire hierarchy of meshes. In this case, we need */
508 /* to set call_refine_interpol_?d to nonzero to make the triple refinement */
509 /* work! */
510 /****************************************************************************/
511 if (mesh->dim == 2) {
512 int i;
513
514 call_refine_interpol_1d = 0;
515
516 for(i = 0; i < mem_info->n_slaves; i++) {
517 MESH *slave = mem_info->slaves[i];
518
519 call_refine_interpol_1d +=
520 count_refine_interpol(
521 slave, AI_get_dof_vec_list(slave), false, &fill_flags);
522 }
523
524 }
525
526 #if DIM_MAX > 2
527 if (mesh->dim == 3) {
528 int i, j;
529
530 call_refine_interpol_1d = 0;
531 call_refine_interpol_2d = 0;
532
533 for(i = 0; i < mem_info->n_slaves; i++) {
534 MESH *slave = mem_info->slaves[i];
535 MESH_MEM_INFO *mem_info_2 = (MESH_MEM_INFO *)slave->mem_info;
536
537 call_refine_interpol_2d +=
538 count_refine_interpol(
539 slave, AI_get_dof_vec_list(slave), false, &fill_flags);
540 if (slave->is_periodic)
541 call_refine_interpol_np_2d +=
542 count_refine_interpol(
543 slave, AI_get_dof_vec_list_np(slave), true, &fill_flags);
544
545 for(j = 0; j < mem_info_2->n_slaves; j++) {
546 MESH *slave_2 = mem_info_2->slaves[j];
547
548 call_refine_interpol_1d +=
549 count_refine_interpol(
550 slave_2, AI_get_dof_vec_list(slave_2), false, &fill_flags);
551 }
552 }
553 }
554 #endif
555 #endif
556 }
557
558 if (mem_info->master) {
559 /****************************************************************************/
560 /* We are on a slave mesh. */
561 /****************************************************************************/
562 int n_slave_elements = mesh->n_elements;
563
564 /* Transfer the refinement marks to the master mesh. */
565 do {
566 mesh_traverse(mesh, 0, FILL_NOTHING|CALL_LEAF_EL, transfer_fct, NULL);
567 mesh_refined = refine(mem_info->master, fill_flags);
568 } while (mesh_refined);
569 return (mesh->n_elements > n_slave_elements) ? MESH_REFINED : 0;
570 #if 0
571 } else {
572 /****************************************************************************/
573 /* We are on a top level master mesh. */
574 /****************************************************************************/
575
576 /* Advance cookies on this mesh and all its slaves. */
577 AI_advance_cookies_rec(mesh);
578 #endif
579 }
580
581 switch(mesh->dim) {
582 case 0:
583 WARNING("No refinement possible for dim == 0!\n");
584 break;
585 case 1:
586 mesh_refined = refine_1d(mesh, fill_flags);
587 break;
588 #if DIM_MAX > 1
589 case 2:
590 mesh_refined = refine_2d(mesh, fill_flags);
591 break;
592 #if DIM_MAX > 2
593 case 3:
594 mesh_refined = refine_3d(mesh, fill_flags);
595 break;
596 #endif
597 #endif
598 default:
599 ERROR_EXIT("Illegal dim during refining!\n");
600 break;
601 }
602
603 if (mesh_refined) {
604 /* Advance cookies on this mesh and all its slaves. */
605 AI_advance_cookies_rec(mesh);
606 }
607
608 return mesh_refined;
609 }
610