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 /* */
7 /* file: submesh_2d.c */
8 /* */
9 /* */
10 /* description: Support for master/slave meshes for master->dim == 2 */
11 /* */
12 /*--------------------------------------------------------------------------*/
13 /* */
14 /* authors: Daniel Koester */
15 /* Institut fuer Mathematik */
16 /* Universitaet Augsburg */
17 /* Universitaetsstr. 14 */
18 /* D-86159 Augsburg, Germany */
19 /* */
20 /* Claus-Justus Heine */
21 /* Abteilung fuer Angewandte Mathematik */
22 /* Albert-Ludwigs-Universitaet Freiburg */
23 /* Hermann-Herder-Str. 10 */
24 /* D-79104 Freiburg im Breisgau, Germany */
25 /* */
26 /* http://www.alberta-fem.de */
27 /* */
28 /* (c) by D. Koester (2005) */
29 /*--------------------------------------------------------------------------*/
30
31
get_slave_elements_rec_2d(MESH * master,MESH * slave,int neigh,EL * m_el,EL * s_el)32 static void get_slave_elements_rec_2d(MESH *master, MESH *slave,
33 int neigh, EL *m_el, EL *s_el)
34 {
35 /* FUNCNAME("get_slave_elements_rec_2d"); */
36 EL_INFO el_info = { NULL, };
37
38
39 if (m_el->child[0]) {
40 if(neigh == 0)
41 get_slave_elements_rec_2d(master, slave, 2, m_el->child[1], s_el);
42 else if (neigh == 1)
43 get_slave_elements_rec_2d(master, slave, 2, m_el->child[0], s_el);
44 else if (neigh == 2) {
45 /* We MUST copy m_el->new_coord here!!! */
46 if (m_el->new_coord) {
47 s_el->new_coord = get_real_d(slave);
48 COPY_DOW(m_el->new_coord, s_el->new_coord);
49 }
50 el_info.mesh = slave;
51 el_info.el = s_el;
52 s_el->mark = 1;
53
54 AI_refine_fct_1d(&el_info, NULL);
55
56 get_slave_elements_rec_2d(master, slave, 0,
57 m_el->child[0], s_el->child[0]);
58 get_slave_elements_rec_2d(master, slave, 1,
59 m_el->child[1], s_el->child[1]);
60 }
61 }
62
63 return;
64 }
65
66
get_slave_elements_2d(MESH * master,MESH * slave,bool (* binding_method)(MESH * master,MACRO_EL * el,int wall,void * data),void * data)67 static void get_slave_elements_2d(MESH *master, MESH *slave,
68 bool (*binding_method)
69 (MESH *master, MACRO_EL *el, int wall,
70 void *data), void *data)
71 {
72 FUNCNAME("get_slave_elements_2d");
73 MACRO_EL *m_mel, *s_mel;
74 int n, i;
75
76 s_mel = slave->macro_els;
77 for(n = 0; n < master->n_macro_el; n++) {
78 m_mel = master->macro_els + n;
79
80 for (i = 0; i < N_NEIGH_2D; i++)
81 if (binding_method(master, m_mel, i, data)) {
82 DEBUG_TEST_EXIT(s_mel,
83 "Ran out of slave macro elements... Wrong meshes?\n");
84
85 get_slave_elements_rec_2d(master, slave,
86 i, m_mel->el, s_mel->el);
87
88 s_mel++;
89 }
90 }
91
92 return;
93 }
94
95
96 /****************************************************************************/
97 /* master interpolation/restriction routines. */
98 /* These call the corresponding refinement/coarsening routines for the slave*/
99 /* mesh patchwise. */
100 /****************************************************************************/
101
master_interpol_2d(DOF_PTR_VEC * m_dpv,RC_LIST_EL * rclist,int mn)102 static void master_interpol_2d(DOF_PTR_VEC *m_dpv, RC_LIST_EL *rclist, int mn)
103 {
104 FUNCNAME("master_interpol_2d");
105 MESH_MEM_INFO *m_mem_info = (MESH_MEM_INFO *)
106 m_dpv->fe_space->admin->mesh->mem_info;
107 int m_n0 = m_dpv->fe_space->admin->n0_dof[EDGE];
108 int m_n = m_dpv->fe_space->admin->mesh->node[EDGE];
109 EL *m_el, *s_el;
110 EL *s_child[2], *m_child[2];
111 int s_n0, s_n, i, j, n_slaves = m_mem_info->n_slaves;
112 MESH *slave = NULL;
113 DOF_PTR_VEC *s_dpv;
114 const NODE_PROJECTION *s_proj;
115
116 /****************************************************************************/
117 /* Retrieve the slave mesh. Rather unelegant, sorry... */
118 /****************************************************************************/
119 for (j = 0; j < n_slaves; j++) {
120 slave = m_mem_info->slaves[j];
121 if(((MESH_MEM_INFO *)slave->mem_info)->slave_binding == m_dpv) break;
122 }
123 DEBUG_TEST_EXIT(j < n_slaves, "Slave mesh not found!\n");
124
125 s_dpv = ((MESH_MEM_INFO *)slave->mem_info)->master_binding;
126 s_n0 = s_dpv->fe_space->admin->n0_dof[CENTER];
127 s_n = slave->node[CENTER];
128
129 /****************************************************************************/
130 /* Run over all patch elements. Check if any edges belong to slave elements.*/
131 /* If so, determine if they should be refined together with the master */
132 /* elements. All necessary DOFs are set. New DOFs are set to NULL, if they */
133 /* should not point to anything. THIS IS NOT JUST A SECURITY MEASURE - IT */
134 /* IS ABSOLUTELY NECESSARY! At the moment, new entries in a DOF vector can */
135 /* not be guaranteed to be clean. */
136 /****************************************************************************/
137
138 for(i = 0; i < mn; i++) {
139 m_el = rclist[i].el_info.el;
140
141 m_child[0] = m_el->child[0];
142 m_child[1] = m_el->child[1];
143
144 /****************************************************************************/
145 /* DOF pointers of the master child elements along the bisection line are */
146 /* set to NULL. */
147 /****************************************************************************/
148 m_dpv->vec[m_child[0]->dof[m_n+1][m_n0]] = NULL;
149 m_dpv->vec[m_child[1]->dof[m_n+0][m_n0]] = NULL;
150
151 for(j = 0; j < N_WALLS_2D; j++) {
152 s_el = (EL *)m_dpv->vec[m_el->dof[m_n+j][m_n0]];
153
154 if(s_el && m_el==(EL *)s_dpv->vec[s_el->dof[s_n][s_n0]]) {
155 EL_INFO el_info = { NULL, };
156
157 if(j == 2) {
158 /****************************************************************************/
159 /* Refine this slave element. Take precautions to preserve the coarse DOFs */
160 /* on master and slave elements, since they will be needed for coarsening. */
161 /****************************************************************************/
162 s_el->mark = MAX(s_el->mark, 1);
163 el_info.el = s_el;
164 el_info.mesh = slave;
165 el_info.master.el = m_el;
166 el_info.master.opp_vertex = j;
167
168 /* Fill active projection on slave rclist for parametric meshes. */
169 s_proj = wall_proj(&rclist[i].el_info, j);
170 if(!s_proj)
171 s_proj = wall_proj(&rclist[i].el_info, -1);
172 el_info.active_projection = s_proj;
173
174 /* DK: Also fill coordinate information on the ad-hoc el_info. This is */
175 /* necessary when we have a non-parametric mesh with projection functions, */
176 /* i.e. we are using the traditional "new_coord" mechanism. */
177
178 COPY_DOW(rclist[i].el_info.coord[0],
179 el_info.coord[0]);
180 COPY_DOW(rclist[i].el_info.coord[1],
181 el_info.coord[1]);
182
183 AI_refine_fct_1d(&el_info, NULL);
184
185 /****************************************************************************/
186 /* Now retrieve the new child elements and set DOF pointers. */
187 /****************************************************************************/
188 s_child[0] = s_el->child[0];
189 s_child[1] = s_el->child[1];
190
191 m_dpv->vec[m_child[0]->dof[m_n+0][m_n0]] = s_child[0];
192 m_dpv->vec[m_child[1]->dof[m_n+1][m_n0]] = s_child[1];
193
194 s_dpv->vec[s_child[0]->dof[s_n][s_n0]] = m_child[0];
195 s_dpv->vec[s_child[1]->dof[s_n][s_n0]] = m_child[1];
196
197 }
198 else if (s_el) {
199 /****************************************************************************/
200 /* If we are not on a refinement edge, set the correct master child */
201 /* pointer to point to the slave element. Edge 2 of child[1-j] should point */
202 /* to the current slave element, and vice versa. */
203 /****************************************************************************/
204 m_dpv->vec[m_child[1-j]->dof[m_n+2][m_n0]] = s_el;
205
206 s_dpv->vec[s_el->dof[s_n][s_n0]] = m_child[1-j];
207 }
208 }
209 /****************************************************************************/
210 /* If there is no slave element on this refinement edge, zero the master */
211 /* child pointers. */
212 /****************************************************************************/
213 else if (!s_el) {
214 if(j == 2) {
215 m_dpv->vec[m_child[0]->dof[m_n+0][m_n0]] = NULL;
216 m_dpv->vec[m_child[1]->dof[m_n+1][m_n0]] = NULL;
217 }
218 else
219 m_dpv->vec[m_child[1-j]->dof[m_n+2][m_n0]] = NULL;
220 }
221 }
222 }
223
224 return;
225 }
226
master_restrict_2d(DOF_PTR_VEC * m_dpv,RC_LIST_EL * rclist,int mn)227 static void master_restrict_2d(DOF_PTR_VEC *m_dpv, RC_LIST_EL *rclist, int mn)
228 {
229 FUNCNAME("master_restrict_2d");
230 MESH_MEM_INFO *m_mem_info = (MESH_MEM_INFO *)
231 m_dpv->fe_space->admin->mesh->mem_info;
232 int m_n0 = m_dpv->fe_space->admin->n0_dof[EDGE];
233 int m_n = m_dpv->fe_space->admin->mesh->node[EDGE];
234 EL *m_el, *s_el, *cm_el;
235 EL *s_child[2], *m_child[2];
236 int s_n0, s_n, i, j, n_slaves = m_mem_info->n_slaves;
237 MESH *slave = NULL;
238 DOF_PTR_VEC *s_dpv;
239
240 /****************************************************************************/
241 /* Retrieve the slave mesh. Rather unelegant, sorry... */
242 /****************************************************************************/
243 for (j = 0; j < n_slaves; j++) {
244 slave = m_mem_info->slaves[j];
245 if(((MESH_MEM_INFO *)slave->mem_info)->slave_binding == m_dpv) break;
246 }
247 DEBUG_TEST_EXIT(j < n_slaves, "Slave mesh not found!\n");
248
249 s_dpv = ((MESH_MEM_INFO *)slave->mem_info)->master_binding;
250 s_n0 = s_dpv->fe_space->admin->n0_dof[CENTER];
251 s_n = slave->node[CENTER];
252
253 /****************************************************************************/
254 /* Run over all patch elements. Check if any edges belong to slave elements.*/
255 /* If so, determine if they should be coarsened together with the master */
256 /* elements. */
257 /****************************************************************************/
258
259 for(i = 0; i < mn; i++) {
260 m_el = rclist[i].el_info.el;
261
262 m_child[0] = m_el->child[0];
263 m_child[1] = m_el->child[1];
264
265 for(j = 0; j < N_WALLS_2D; j++) {
266 s_el = (EL *)m_dpv->vec[m_el->dof[m_n+j][m_n0]];
267
268 if(s_el) {
269 cm_el = (EL *)s_dpv->vec[s_el->dof[s_n][s_n0]];
270
271 if (cm_el==m_child[0] || cm_el==m_child[1] || cm_el==m_el) {
272 /****************************************************************************/
273 /* Set the slave pointer to point to the parent master element. */
274 /****************************************************************************/
275 s_dpv->vec[s_el->dof[s_n][s_n0]] = m_el;
276
277 if (j==2) {
278 /****************************************************************************/
279 /* Coarsen this slave element. Old parent DOFs still exist due to the */
280 /* preserve_coarse_dofs-setting during refinement. */
281 /****************************************************************************/
282 EL_INFO el_info = { NULL, };
283
284 s_child[0] = s_el->child[0];
285 s_child[1] = s_el->child[1];
286
287 s_child[0]->mark = -1;
288 s_child[1]->mark = -1;
289
290 el_info.mesh = slave;
291 el_info.el = s_el;
292 AI_coarse_fct_1d(&el_info, NULL);
293 }
294 }
295 }
296 }
297 }
298
299 return;
300 }
301
302
join_elements_recursive_2d(const MESH * master,const MESH * slave,const DOF_ADMIN * m_admin,const DOF_ADMIN * s_admin,const DOF_PTR_VEC * m_dpv,const DOF_PTR_VEC * s_dpv,const int subsimplex,const EL * m_el,const EL * s_el)303 static void join_elements_recursive_2d(const MESH *master,
304 const MESH *slave,
305 const DOF_ADMIN *m_admin,
306 const DOF_ADMIN *s_admin,
307 const DOF_PTR_VEC *m_dpv,
308 const DOF_PTR_VEC *s_dpv,
309 const int subsimplex,
310 const EL *m_el,
311 const EL *s_el)
312 {
313 FUNCNAME("join_elements_recursive_2d");
314
315 s_dpv->vec[s_el->dof[slave->node[CENTER]]
316 [s_admin->n0_dof[CENTER]]] = (void *)m_el;
317
318 m_dpv->vec[m_el->dof[master->node[EDGE] + subsimplex]
319 [m_admin->n0_dof[EDGE]]] = (void *)s_el;
320
321 if (m_el->child[0]) {
322 if (subsimplex == 2) { /* this is the refinement edge */
323 DEBUG_TEST_EXIT(s_el->child[0],
324 "Could not find slave children!\n");
325
326 join_elements_recursive_2d(master, slave, m_admin, s_admin,
327 m_dpv, s_dpv,
328 0, m_el->child[0], s_el->child[0]);
329 join_elements_recursive_2d(master, slave, m_admin, s_admin,
330 m_dpv, s_dpv,
331 1, m_el->child[1], s_el->child[1]);
332 }
333 else
334 join_elements_recursive_2d(master, slave, m_admin, s_admin,
335 m_dpv, s_dpv,
336 2, m_el->child[1 - subsimplex], s_el);
337 }
338
339 return;
340 }
341
342 /* If we have a periodic mesh then we need to transfer the wall
343 * transformations to the slave mesh. Of course, the slave can only
344 * inherit transformations which make sense, i.e. do not map vertices
345 * away from the child.
346 *
347 * Wall transformations must be set before calling
348 * compute_neigh_fast()
349 *
350 * We may produce duplicates here, but that does not matter, its just
351 * a slight performance hit.
352 */
transfer_wall_trafos_2d(MESH * master,MACRO_DATA * sdata,int * vert_ind)353 static void transfer_wall_trafos_2d(MESH *master,
354 MACRO_DATA *sdata,
355 int *vert_ind)
356 {
357 int (*master_wall_vtx_trafos)[N_VERTICES(DIM_MAX-1)][2];
358 int (*slave_wall_vtx_trafos)[N_VERTICES(DIM_MAX-1)][2] = NULL;
359 int nwt;
360 int nswt, v, i;
361
362 nwt = _AI_compute_macro_wall_trafos(master, &master_wall_vtx_trafos);
363
364 nswt = 0;
365 for (i = 0; i < nwt; i++) {
366 for (v = 0; v < N_VERTICES_1D; v++) {
367 if (vert_ind[master_wall_vtx_trafos[i][v][0]] >= 0 &&
368 vert_ind[master_wall_vtx_trafos[i][v][1]] >= 0) {
369 ++nswt;
370 }
371 }
372 }
373
374 if (nswt > 0) {
375 slave_wall_vtx_trafos = (int (*)[N_VERTICES(DIM_MAX-1)][2])
376 MEM_ALLOC(nswt*sizeof(*slave_wall_vtx_trafos), char);
377
378 nswt = 0;
379 for (i = 0; i < nwt; i++) {
380 for (v = 0; v < N_VERTICES_1D; v++) {
381 int from = vert_ind[master_wall_vtx_trafos[i][v][0]];
382 int to = vert_ind[master_wall_vtx_trafos[i][v][1]];
383
384 if (from >= 0 && to >= 0) {
385 slave_wall_vtx_trafos[nswt][0][0] = from;
386 slave_wall_vtx_trafos[nswt][0][1] = to;
387 ++nswt;
388 }
389 }
390 }
391 sdata->n_wall_vtx_trafos = nswt;
392 sdata->wall_vtx_trafos = slave_wall_vtx_trafos;
393
394 sdata->el_wall_vtx_trafos =
395 MEM_ALLOC(sdata->n_macro_elements*N_WALLS_1D, int);
396 _AI_compute_element_wall_transformations(sdata);
397 }
398
399 MEM_FREE(master_wall_vtx_trafos, nwt*sizeof(*master_wall_vtx_trafos), char);
400 }
401
402 /****************************************************************************/
403 /* get_submesh_2d(master, name, binding_method): returns a 1D submesh of */
404 /* master based on the information given by the user routine */
405 /* binding_method(). */
406 /****************************************************************************/
407
408 static MESH *
get_submesh_2d(MESH * master,const char * name,bool (* binding_method)(MESH * master,MACRO_EL * el,int edge,void * data),void * data)409 get_submesh_2d(MESH *master, const char *name,
410 bool (*binding_method)(MESH *master,
411 MACRO_EL *el, int edge,
412 void *data),
413 void *data)
414 {
415 FUNCNAME("get_submesh_2d");
416 MACRO_DATA s_data = { 0 /* dim */, };
417 MESH_MEM_INFO *m_info, *s_info;
418 MESH *slave = NULL;
419 const FE_SPACE *slave_space, *master_space;
420 DOF_PTR_VEC *slave_to_master_binding;
421 DOF_PTR_VEC *master_to_slave_binding;
422 int s_n_dof[N_NODE_TYPES] = { 0, };
423 int m_n_dof[N_NODE_TYPES] = { 0, };
424 MACRO_EL *m_mel, *s_mel;
425 const DOF_ADMIN *m_admin, *s_admin;
426 int i, j, k, n, ne = 0, nv = 0, *vert_ind = NULL, index;
427 char new_name[1024];
428
429 m_info = ((MESH_MEM_INFO *)master->mem_info);
430
431 /****************************************************************************/
432 /* Count all needed vertices and elements. */
433 /****************************************************************************/
434
435 s_data.dim = 1;
436 s_data.coords = MEM_ALLOC(master->n_vertices, REAL_D); /* resized later! */
437
438 vert_ind = MEM_ALLOC(master->n_vertices, int);
439 for (i = 0; i < master->n_vertices; i++)
440 vert_ind[i] = -1;
441
442 for (n = 0; n < master->n_macro_el; n++) {
443 m_mel = master->macro_els + n;
444
445 for (i = 0; i < N_WALLS_2D; i++)
446 if (binding_method(master, m_mel, i, data)) {
447 ne++;
448
449 for (j = 0; j < N_VERTICES_2D; j++)
450 if (j != i) {
451 /* Make use of the slave->mem_info->coords vector to get vertex indices. */
452 index = m_mel->coord[j] - &m_info->coords[0];
453
454 if (vert_ind[index] < 0) {
455 vert_ind[index] = nv;
456 for(k = 0; k < DIM_OF_WORLD; k++)
457 s_data.coords[nv][k] = m_info->coords[index][k];
458
459 nv++;
460 }
461 }
462 }
463 }
464
465 /****************************************************************************/
466 /* Allocate the needed amount of macro elements and coordinates. */
467 /* Fill element and coordinate information. */
468 /****************************************************************************/
469
470 TEST_EXIT(nv,"Bad mesh: no vertices counted!\n");
471 TEST_EXIT(ne,"Bad mesh: no elements counted!\n");
472
473 s_data.n_total_vertices = nv;
474 s_data.n_macro_elements = ne;
475 s_data.coords = MEM_REALLOC(s_data.coords, master->n_vertices,
476 nv, REAL_D);
477 s_data.mel_vertices = MEM_ALLOC(ne * N_NEIGH_1D, int);
478 ne = 0;
479
480 for (n = 0; n < master->n_macro_el; n++) {
481 m_mel = master->macro_els + n;
482
483 for (i = 0; i < N_WALLS_2D; i++)
484 if (binding_method(master, m_mel, i, data)) {
485 for (j = 0; j < N_VERTICES_2D; j++)
486 if (j != i) {
487 index = m_mel->coord[j] - &m_info->coords[0];
488 nv = vert_ind[index];
489
490 s_data.mel_vertices[N_VERTICES_1D * ne + (j + 2 - i) % 3] = nv;
491 }
492
493 ne++;
494 }
495 }
496
497 /* Wall transformations must be set before calling
498 * compute_neigh_fast(). Slaves of periodic meshes need not be
499 * periodic.
500 */
501 if (master->is_periodic) {
502 transfer_wall_trafos_2d(master, &s_data, vert_ind);
503 }
504
505 compute_neigh_fast(&s_data);
506
507 /* Assign all boundaries a no-boundary type (also allocates .boundary) */
508 default_boundary(&s_data, INTERIOR, true);
509
510 /* Now correct the boundary values if the master vertex has nonzero type. */
511 ne = 0;
512 for (ne = n = 0; n < master->n_macro_el; n++) {
513 m_mel = master->macro_els + n;
514
515 for (i = 0; i < N_WALLS_2D; i++)
516 if (binding_method(master, m_mel, i, data)) {
517 for (j = 0; j < N_NEIGH_1D; j++) {
518 if (s_data.neigh[N_NEIGH_1D * ne + j] < 0) {
519 BNDRY_FLAGS mask;
520 BNDRY_FLAGS_INIT(mask);
521 BNDRY_FLAGS_SET(mask, m_mel->wall_bound[i]);
522 BNDRY_FLAGS_XOR(mask, m_mel->vertex_bound[(i+2-j) % N_VERTICES_2D]);
523 int bound = BNDRY_FLAGS_FFBB(mask);
524 if (bound > INTERIOR &&
525 (s_data.boundary[N_NEIGH_1D * ne + j] == INTERIOR ||
526 (unsigned)bound < s_data.boundary[N_NEIGH_1D * ne + j])) {
527 s_data.boundary[N_NEIGH_1D * ne + j] = bound;
528 }
529 }
530 }
531 ne++;
532 }
533 }
534
535 /* Assign all remaining boundaries a boundary type of 1. Should we
536 * make this configurable? Or use (N_BNDRY_BITS-1) in order not to
537 * interfere with existing boundary classifications?
538 */
539 default_boundary(&s_data, DIRICHLET, false);
540
541 /****************************************************************************/
542 /* Allocate a submesh. */
543 /****************************************************************************/
544
545 if(!name) {
546 static int count_2d = 1;
547
548 sprintf(new_name, "Submesh %d of %s", count_2d, master->name);
549 name = new_name;
550
551 count_2d++;
552 }
553
554 slave = GET_MESH(1, name, &s_data, NULL, NULL);
555
556 /****************************************************************************/
557 /* Clean up. */
558 /****************************************************************************/
559
560 nv = s_data.n_total_vertices;
561 ne = s_data.n_macro_elements;
562
563 MEM_FREE(s_data.coords, nv, REAL_D);
564 MEM_FREE(s_data.mel_vertices, ne*N_VERTICES_1D, int);
565 MEM_FREE(s_data.neigh, ne*N_NEIGH_1D, int);
566 MEM_FREE(s_data.opp_vertex, ne*N_NEIGH_1D, int);
567 MEM_FREE(s_data.boundary, ne*N_NEIGH_1D, S_CHAR);
568 MEM_FREE(vert_ind, master->n_vertices, int);
569
570 /****************************************************************************/
571 /* Fill more slave elements, if the master mesh was already refined. */
572 /****************************************************************************/
573
574 get_slave_elements_2d(master, slave, binding_method, data);
575
576 /****************************************************************************/
577 /* Allocate special FE spaces for the slave. */
578 /****************************************************************************/
579
580 s_n_dof[CENTER] = 1;
581
582 slave_space = get_dof_space(slave, "Center dof fe_space", s_n_dof,
583 ADM_PRESERVE_COARSE_DOFS);
584
585 slave_to_master_binding = get_dof_ptr_vec("Slave - master pointers",
586 slave_space);
587
588 /****************************************************************************/
589 /* Allocate special FE spaces for master. */
590 /****************************************************************************/
591
592 m_n_dof[EDGE] = 1;
593 master_space = get_dof_space(master, "Edge dof fe_space", m_n_dof,
594 ADM_PRESERVE_COARSE_DOFS);
595
596 #if ALBERTA_DEBUG == 1
597 check_mesh(slave);
598 #endif
599
600 /****************************************************************************/
601 /* Allocate special DOF_PTR_VECs for both master and slave. These serve to */
602 /* help find the corresponding subsimplex to each boundary master simplex */
603 /* during refinement and vice versa. */
604 /****************************************************************************/
605
606 master_to_slave_binding = get_dof_ptr_vec("Master - slave pointers",
607 master_space);
608
609 master_to_slave_binding->refine_interpol = master_interpol_2d;
610 master_to_slave_binding->coarse_restrict = master_restrict_2d;
611
612 /****************************************************************************/
613 /* Set the special pointers in the MESH_MEM_INFO components of both master */
614 /* and slave grids. */
615 /****************************************************************************/
616
617 s_info = (MESH_MEM_INFO *)slave->mem_info;
618 s_info->master = master;
619 s_info->slave_binding = master_to_slave_binding;
620 s_info->master_binding = slave_to_master_binding;
621
622 m_info->slaves = MEM_REALLOC(m_info->slaves,
623 m_info->n_slaves,
624 m_info->n_slaves + 1,
625 MESH *);
626 m_info->slaves[m_info->n_slaves] = slave;
627
628 m_info->n_slaves++;
629
630 /****************************************************************************/
631 /* Set the element pointer vec entries to the correct values. */
632 /* This assumes that slave macro elements were allocated in the order given */
633 /* by the loop below. */
634 /****************************************************************************/
635
636 m_admin = master_to_slave_binding->fe_space->admin;
637 s_admin = slave_to_master_binding->fe_space->admin;
638
639 FOR_ALL_DOFS(s_admin, slave_to_master_binding->vec[dof] = NULL);
640 FOR_ALL_DOFS(m_admin, master_to_slave_binding->vec[dof] = NULL);
641
642 s_mel = slave->macro_els;
643 for(n = 0; n < master->n_macro_el; n++) {
644 int ov;
645
646 m_mel = master->macro_els + n;
647
648 for (ov = 0; ov < N_NEIGH_2D; ov++)
649 if (binding_method(master, m_mel, ov, data)) {
650 int v0, v1;
651
652 DEBUG_TEST_EXIT(s_mel,
653 "Ran out of slave macro elements... Wrong meshes?\n");
654
655 /* Here we take care of node projection function transfer. */
656 if(m_mel->projection[ov+1])
657 s_mel->projection[0] = m_mel->projection[ov+1];
658 else
659 s_mel->projection[0] = m_mel->projection[0];
660
661 join_elements_recursive_2d(master, slave, m_admin, s_admin,
662 master_to_slave_binding,
663 slave_to_master_binding,
664 ov, m_mel->el, s_mel->el);
665
666 s_mel->master.macro_el = m_mel;
667 s_mel->master.opp_vertex = ov;
668
669 /* also store the boundary classification of the master
670 * vertices and edges.
671 */
672 v0 = vertex_of_wall_2d[ov][0];
673 v1 = vertex_of_wall_2d[ov][1];
674 if (m_info->master) {
675 /* slave hierarchy. */
676 BNDRY_FLAGS_CPY(s_mel->master.vertex_bound[0],
677 m_mel->master.vertex_bound[v0]);
678 BNDRY_FLAGS_CPY(s_mel->master.np_vertex_bound[0],
679 m_mel->master.np_vertex_bound[v0]);
680 BNDRY_FLAGS_CPY(s_mel->master.vertex_bound[1],
681 m_mel->master.vertex_bound[v1]);
682 BNDRY_FLAGS_CPY(s_mel->master.np_vertex_bound[1],
683 m_mel->master.np_vertex_bound[v1]);
684 BNDRY_FLAGS_CPY(s_mel->master.edge_bound[0],
685 m_mel->master.edge_bound[ov]);
686 BNDRY_FLAGS_CPY(s_mel->master.np_edge_bound[0],
687 m_mel->master.np_edge_bound[ov]);
688 } else {
689 BNDRY_FLAGS_CPY(s_mel->master.vertex_bound[0],
690 m_mel->vertex_bound[v0]);
691 BNDRY_FLAGS_CPY(s_mel->master.np_vertex_bound[0],
692 m_mel->np_vertex_bound[v0]);
693 BNDRY_FLAGS_CPY(s_mel->master.vertex_bound[1],
694 m_mel->vertex_bound[v1]);
695 BNDRY_FLAGS_CPY(s_mel->master.np_vertex_bound[1],
696 m_mel->np_vertex_bound[v1]);
697 BNDRY_FLAGS_INIT(s_mel->master.edge_bound[0]);
698 BNDRY_FLAGS_INIT(s_mel->master.np_edge_bound[0]);
699 if (m_mel->neigh_vertices[ov][0] == -1) {
700 BNDRY_FLAGS_SET(s_mel->master.edge_bound[0],
701 m_mel->wall_bound[ov]);
702 }
703 BNDRY_FLAGS_SET(s_mel->master.np_edge_bound[0],
704 m_mel->wall_bound[ov]);
705 }
706
707 s_mel++;
708 }
709 }
710
711 return slave;
712 }
713
714