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.c */
8 /* */
9 /* */
10 /* description: Support for master/slave meshes */
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 /* Universitaet Freiburg */
23 /* Hermann-Herder-Strasse 10 */
24 /* 79104 Freiburg, Germany */
25 /* */
26 /* http://www.alberta-fem.de */
27 /* */
28 /* (c) by D. Koester (2004), */
29 /* C.-J. Heine (2006-2009). */
30 /*--------------------------------------------------------------------------*/
31
32 #include "alberta_intern.h"
33 #include "alberta.h"
34
35 #include "submesh_1d.c"
36 #if DIM_MAX > 1
37 #include "submesh_2d.c"
38 #if DIM_MAX > 2
39 #include "submesh_3d.c"
40 #endif
41 #endif
42
43 static inline
bind_one_mel(MESH * master,MACRO_EL * m_mel,int wall,DOF_PTR_VEC * s_dpv,MESH * slave,MACRO_EL * s_mel,DOF_PTR_VEC * m_dpv)44 void bind_one_mel(MESH *master, MACRO_EL *m_mel, int wall, DOF_PTR_VEC *s_dpv,
45 MESH *slave, MACRO_EL *s_mel, DOF_PTR_VEC *m_dpv)
46 {
47 int m_dim = master->dim;
48
49 /* Here we take care of node projection function transfer. */
50 if (m_dim > 1) {
51 if (m_mel->projection[wall+1])
52 s_mel->projection[0] = m_mel->projection[wall+1];
53 else
54 s_mel->projection[0] = m_mel->projection[0];
55 }
56
57 if (m_dim == 1)
58 join_elements_recursive_1d(master, slave,
59 m_dpv->fe_space->admin,
60 s_dpv->fe_space->admin,
61 m_dpv, s_dpv,
62 wall, m_mel->el, s_mel->el);
63 #if DIM_MAX > 1
64 else if (m_dim == 2)
65 join_elements_recursive_2d(master, slave,
66 m_dpv->fe_space->admin,
67 s_dpv->fe_space->admin,
68 m_dpv, s_dpv,
69 wall, m_mel->el, s_mel->el);
70 #endif
71 #if DIM_MAX >= 3
72 else if (m_dim == 3)
73 join_elements_recursive_3d(master, slave,
74 m_dpv->fe_space->admin,
75 s_dpv->fe_space->admin,
76 m_dpv, s_dpv,
77 wall, m_mel->el, s_mel->el,
78 m_mel->orientation, m_mel->el_type);
79 #endif
80 }
81
82 /* Bind a given mesh (e.g. read from disk by read_mesh) to a master
83 * mesh. SLAVE and MASTER should be compatible, otherwise the routine
84 * will catch an error.
85 */
bind_submesh(MESH * master,MESH * slave,bool (* binding_method)(MESH * master,MACRO_EL * el,int wall,void * data),void * data)86 void bind_submesh(MESH *master,
87 MESH *slave,
88 bool (*binding_method)(MESH *master, MACRO_EL *el,
89 int wall, void *data),
90 void *data)
91 {
92 FUNCNAME("bind_submesh");
93 MESH_MEM_INFO *s_info, *m_info;
94 MACRO_EL *m_mel, *s_mel;
95 const DOF_ADMIN *admin = NULL;
96 const FE_SPACE *m_space, *s_space;
97 DOF_PTR_VEC *m_dpv, *s_dpv;
98 int n_dof[N_NODE_TYPES] = { 0, };
99 int i, j, n, m_dim;
100
101 /****************************************************************************/
102 /* Do the checks for obvious user errors. */
103 /****************************************************************************/
104
105 TEST_EXIT(master, "No master mesh given!\n");
106 TEST_EXIT(master->dim > 0, "Master mesh has dim == 0!\n");
107 m_dim = master->dim;
108
109 /****************************************************************************/
110 /* Set the mem_info components in master and slave. */
111 /****************************************************************************/
112
113 s_info = (MESH_MEM_INFO *)slave->mem_info;
114 m_info = (MESH_MEM_INFO *)master->mem_info;
115
116 m_info->slaves = MEM_REALLOC(m_info->slaves,
117 m_info->n_slaves,
118 m_info->n_slaves + 1,
119 MESH *);
120
121 m_info->slaves[m_info->n_slaves] = slave;
122
123 m_info->n_slaves++;
124
125 slave->trace_id = m_info->next_trace_id++;
126
127 s_info->master = master;
128
129 /* Check for the correct DOF_ADMINs for slave and master. */
130
131 n_dof[CENTER] = 1;
132
133 for (i = 0; i < slave->n_dof_admin; i++) {
134 admin = slave->dof_admin[i];
135 for (j = 0; j < N_NODE_TYPES; j++) {
136 if (admin->n_dof[j] != n_dof[j]) goto bad_admin;
137 }
138 if (!(admin->flags == ADM_PRESERVE_COARSE_DOFS)) goto bad_admin;
139 break;
140 bad_admin:
141 admin = NULL;
142 }
143 TEST_EXIT(admin, "Slave mesh does not seem to have had a master!\n");
144 s_space = get_dof_space(slave, "Center FE_SPACE", n_dof,
145 ADM_PRESERVE_COARSE_DOFS);
146
147 n_dof[CENTER] = 0;
148 switch(m_dim) {
149 case 1:
150 n_dof[VERTEX] = 1;
151 break;
152 case 2:
153 n_dof[EDGE] = 1;
154 break;
155 case 3:
156 n_dof[FACE] = 1;
157 break;
158 }
159
160 for (i = 0; i < master->n_dof_admin; i++) {
161 admin = master->dof_admin[i];
162 for (j = 0; j < N_NODE_TYPES; j++) {
163 if (admin->n_dof[j] != n_dof[j]) goto bad_admin2;
164 }
165 if (!(admin->flags == ADM_PRESERVE_COARSE_DOFS)) goto bad_admin2;
166 break;
167 bad_admin2:
168 admin = NULL;
169 }
170 TEST_EXIT(admin,"Given master mesh does not seem to have had slaves!\n");
171 m_space = get_dof_space(master, "Wall FE_SPACE", n_dof,
172 ADM_PRESERVE_COARSE_DOFS);
173
174 /* Allocate element pointer vectors. */
175
176 s_dpv = get_dof_ptr_vec("Slave - master pointers", s_space);
177 s_info->master_binding = s_dpv;
178
179 m_dpv = get_dof_ptr_vec("Master - slave pointers", m_space);
180 s_info->slave_binding = m_dpv;
181
182 switch(m_dim) {
183 case 1:
184 m_dpv->refine_interpol = master_interpol_1d;
185 m_dpv->coarse_restrict = master_restrict_1d;
186 break;
187 #if DIM_MAX > 1
188 case 2:
189 m_dpv->refine_interpol = master_interpol_2d;
190 m_dpv->coarse_restrict = master_restrict_2d;
191 break;
192 #if DIM_MAX >= 3
193 case 3:
194 m_dpv->refine_interpol = master_interpol_3d;
195 m_dpv->coarse_restrict = master_restrict_3d;
196 #endif
197 #endif
198 }
199
200 /****************************************************************************/
201 /* Set the element pointer vec entries to the correct values. */
202 /* This assumes that slave macro elements were allocated in the order given */
203 /* by the loop below. */
204 /****************************************************************************/
205
206 FOR_ALL_DOFS(s_dpv->fe_space->admin, s_dpv->vec[dof] = NULL);
207 FOR_ALL_DOFS(m_dpv->fe_space->admin, m_dpv->vec[dof] = NULL);
208
209 if (binding_method == NULL) {
210 /* We assume here that the two meshes are already chained on the
211 * macro element level.
212 */
213 for (n = 0; n < slave->n_macro_el; n++) {
214 s_mel = slave->macro_els + n;
215 m_mel = s_mel->master.macro_el;
216 i = s_mel->master.opp_vertex;
217 TEST_EXIT(m_mel != NULL,
218 "Meshes are not chained on the macro-element level.\n");
219 TEST_EXIT(i >= 0,
220 "Garbled slave->master binding (macro-element level).\n");
221
222 bind_one_mel(master, m_mel, i, s_dpv, slave, s_mel, m_dpv);
223 }
224 } else {
225 MACRO_EL *s_mel = slave->macro_els;
226 MACRO_EL *s_mel_end = s_mel + slave->n_macro_el;
227 for (n = 0; n < master->n_macro_el; n++) {
228 m_mel = master->macro_els + n;
229
230 for (i = 0; i < N_NEIGH(m_dim); i++)
231 if (binding_method(master, m_mel, i, data)) {
232 TEST_EXIT(s_mel < s_mel_end,
233 "Ran out of slave macro elements... Wrong meshes?\n");
234
235 bind_one_mel(master, m_mel, i, s_dpv, slave, s_mel, m_dpv);
236
237 /* Set links for FILL_MASTER_INFO during mesh-traversal */
238 s_mel->master.macro_el = m_mel;
239 s_mel->master.opp_vertex = i;
240
241 s_mel++;
242 }
243 }
244 }
245
246 free_fe_space(s_space);
247 free_fe_space(m_space);
248 }
249
250 /****************************************************************************/
251 /* read_submesh_gen(read_xdr, master, slave_filename, binding_method, */
252 /* init_boundary): */
253 /* Read a slave mesh from file "slave_name" and bind it to "master". Assumes*/
254 /* that master and slave were written at the same time. Available for native*/
255 /* and XDR formats. Return a pointer to the slave mesh. */
256 /****************************************************************************/
257
258 static
read_submesh_gen(int read_xdr,MESH * master,const char * slave_filename,bool (* binding_method)(MESH * master,MACRO_EL * el,int wall,void * data),void * data)259 MESH *read_submesh_gen(int read_xdr, MESH *master,
260 const char *slave_filename,
261 bool (*binding_method)(MESH *master, MACRO_EL *el,
262 int wall, void *data),
263 void *data)
264 {
265 FUNCNAME("read_submesh_gen");
266 MESH *slave = NULL;
267
268 /****************************************************************************/
269 /* Do the checks for obvious user errors. */
270 /****************************************************************************/
271
272 TEST_EXIT(master, "No master mesh given!\n");
273 TEST_EXIT(master->dim > 0, "Master mesh has dim == 0!\n");
274 TEST_EXIT(slave_filename, "No filename for the slave mesh given!\n");
275 TEST_EXIT(binding_method, "No binding method given!\n");
276
277 /****************************************************************************/
278 /* Read the mesh from the file. Do not use the time value pointer. */
279 /****************************************************************************/
280 if (read_xdr)
281 slave = read_mesh_xdr(slave_filename, NULL, NULL, NULL);
282 else
283 slave = read_mesh(slave_filename, NULL, NULL, NULL);
284
285 bind_submesh(master, slave, binding_method, data);
286
287 return slave;
288 }
289
290 /****************************************************************************/
291 /* Interface to the maintainers */
292 /****************************************************************************/
293
294 /****************************************************************************/
295 /* AI_check_slavery(master): Do a consistency check for all submeshes. */
296 /****************************************************************************/
297
AI_check_slavery(MESH * master)298 void AI_check_slavery(MESH *master)
299 {
300 FUNCNAME("AI_check_slavery");
301 MESH_MEM_INFO *m_mem_info, *s_mem_info;
302 int i, k, n_slaves, slave_el_count;
303 MESH *slave;
304 DOF_PTR_VEC *m_dpv, *s_dpv;
305 TRAVERSE_STACK *stack;
306 const DOF_ADMIN *m_admin, *s_admin;
307 const EL_INFO *m_el_info, *s_el_info;
308 const EL *m_el, *s_el;
309
310 if (!master) {
311 MSG("No mesh provided!\n");
312 return;
313 }
314 TEST_EXIT(m_mem_info = (MESH_MEM_INFO *)(master->mem_info),
315 "No memory management present for \"%s\"!\n", master->name);
316 n_slaves = m_mem_info->n_slaves;
317
318 if (!n_slaves) {
319 INFO(4,4,"Mesh \"%d\" has no slaves.\n", master->name);
320 return;
321 }
322
323 stack = get_traverse_stack();
324
325 /****************************************************************************/
326 /* Run over all slave meshes. */
327 /****************************************************************************/
328 for (k = 0; k < n_slaves; k++) {
329 slave = m_mem_info->slaves[k];
330 TEST_EXIT(slave,"Slave mesh no. %d not found!\n", k);
331 INFO(6,6,"Analysing slave \"%s\"...\n", slave->name);
332
333 TEST_EXIT(slave->dim + 1 == master->dim,
334 "Bad dimension of slave!\n");
335
336 TEST_EXIT(s_mem_info = (MESH_MEM_INFO *)(slave->mem_info),
337 "No memory management present for slave!\n");
338
339 TEST_EXIT(s_mem_info->master == master,
340 "Wrong mem_info->master pointer on slave!\n");
341 TEST_EXIT(m_dpv = s_mem_info->slave_binding,
342 "No binding vector to slave present!\n");
343 TEST_EXIT(s_dpv = s_mem_info->master_binding,
344 "No binding vector to master present!\n");
345
346 INFO(8,8,"Slave mesh has %d subslaves.\n", s_mem_info->n_slaves);
347
348 m_admin = m_dpv->fe_space->admin;
349 s_admin = s_dpv->fe_space->admin;
350
351 INFO(10,10,"Current master leaf elements:\n");
352 m_el_info = traverse_first(stack, master, -1, CALL_LEAF_EL);
353 while(m_el_info) {
354 INFO(10,10,"%d\n", INDEX(m_el_info->el));
355 m_el_info = traverse_next(stack, m_el_info);
356 }
357
358 INFO(10,10,"Current slave leaf elements:\n");
359 s_el_info = traverse_first(stack, slave, -1, CALL_LEAF_EL);
360 while(s_el_info) {
361 INFO(10,10,"%d\n", INDEX(s_el_info->el));
362 s_el_info = traverse_next(stack, s_el_info);
363 }
364
365 /****************************************************************************/
366 /* Run over the slave mesh and check correspondance to master. */
367 /****************************************************************************/
368 slave_el_count = 0;
369 s_el_info = traverse_first(stack, slave, -1, CALL_EVERY_EL_PREORDER);
370 while(s_el_info) {
371 slave_el_count++;
372 s_el = s_el_info->el;
373 INFO(10,10,"Analysing slave el %d...\n", INDEX(s_el));
374 if (!IS_LEAF_EL(s_el))
375 INFO(10,10,"(Child elements: %d, %d)\n", INDEX(s_el->child[0]),
376 INDEX(s_el->child[1]));
377 TEST_EXIT(m_el = (EL *) s_dpv->vec[s_el->dof[slave->node[CENTER]]
378 [s_admin->n0_dof[CENTER]]],
379 "Slave element %d does not point to a master element!\n", INDEX(s_el));
380 INFO(10,10,"slave el %d points to master el %d\n",
381 INDEX(s_el), INDEX(m_el));
382
383 for (i = 0; i < N_NEIGH(master->dim); i++) {
384 if (master->dim == 2) {
385 if (s_el == (EL *)m_dpv->vec[m_el->dof[master->node[EDGE]+i]
386 [m_admin->n0_dof[EDGE]]])
387 break;
388 }
389 else
390 if (s_el == (EL *)m_dpv->vec[m_el->dof[master->node[FACE]+i]
391 [m_admin->n0_dof[FACE]]])
392 break;
393 }
394 TEST_EXIT(i < N_NEIGH(master->dim),
395 "Master element %d does not point back to slave element %d!\n",
396 INDEX(m_el), INDEX(s_el));
397
398 s_el_info = traverse_next(stack, s_el_info);
399 }
400 if (slave_el_count < slave->n_hier_elements)
401 ERROR_EXIT("slave element count == %d < %d == slave->n_elements!\n",
402 slave_el_count, slave->n_elements);
403 if (slave_el_count > slave->n_hier_elements)
404 ERROR_EXIT("slave element count == %d > %d == slave->n_elements!\n",
405 slave_el_count, slave->n_elements);
406
407 /****************************************************************************/
408 /* Run over the current master mesh and check correspondance to slave. */
409 /****************************************************************************/
410 m_el_info = traverse_first(stack, master, -1,
411 CALL_EVERY_EL_PREORDER|FILL_ORIENTATION);
412 while(m_el_info) {
413 m_el = m_el_info->el;
414 INFO(10,10,"Analysing master el %d...\n", INDEX(m_el));
415 if (!IS_LEAF_EL(m_el))
416 INFO(10,10,"(Child elements: %d, %d)\n", INDEX(m_el->child[0]),
417 INDEX(m_el->child[1]));
418
419 for (i = 0; i < N_NEIGH(master->dim); i++) {
420 if (master->dim == 2)
421 s_el = (EL *)m_dpv->vec[m_el->dof[master->node[EDGE]+i]
422 [m_admin->n0_dof[EDGE]]];
423 else
424 s_el = (EL *)m_dpv->vec[m_el->dof[master->node[FACE]+i]
425 [m_admin->n0_dof[FACE]]];
426 if (s_el) {
427 INFO(10,10,"master el %d, subsimplex %d, points to slave el %d\n",
428 INDEX(m_el), i, INDEX(s_el));
429
430 if (IS_LEAF_EL(m_el)) {
431 TEST_EXIT(m_el == (EL *)s_dpv->vec[s_el->dof[slave->node[CENTER]]
432 [s_admin->n0_dof[CENTER]]],
433 "Slave element %d does not point back to master element %d!\n",
434 INDEX(s_el), INDEX(m_el));
435 }
436 }
437 }
438
439 m_el_info = traverse_next(stack, m_el_info);
440 }
441 }
442
443 INFO(4,4,"No errors found.\n");
444 free_traverse_stack(stack);
445
446 return;
447 }
448
449 /****************************************************************************/
450 /* Interface to the user */
451 /****************************************************************************/
452
453
454 /****************************************************************************/
455 /* get_submesh(master, name, init_leaf_data, binding_method): */
456 /* The main allocation routine for getting slave meshes. */
457 /****************************************************************************/
458
get_submesh(MESH * master,const char * name,bool (* binding_method)(MESH * master,MACRO_EL * el,int wall,void * data),void * data)459 MESH *get_submesh(MESH *master, const char *name,
460 bool (*binding_method)(MESH *master,
461 MACRO_EL *el, int wall,
462 void *data),
463 void *data)
464 {
465 FUNCNAME("get_submesh");
466 MESH *slave = NULL;
467
468 /****************************************************************************/
469 /* Do the checks for obvious user errors. */
470 /****************************************************************************/
471
472 TEST_EXIT(master,"No master mesh specified!\n");
473
474 TEST_EXIT(master->dim > 0, "Does not make sense for dim 0 master meshes!\n");
475
476 TEST_EXIT(binding_method, "Parameter 'binding_method' must be nonzero!\n");
477
478 /****************************************************************************/
479 /* Do the dimension dependent stuff for the slave. */
480 /****************************************************************************/
481
482 if (master->dim == 1)
483 slave = get_submesh_1d(master, name, binding_method, data);
484 #if DIM_MAX > 1
485 else if (master->dim == 2)
486 slave = get_submesh_2d(master, name, binding_method, data);
487 #if DIM_MAX > 2
488 else
489 slave = get_submesh_3d(master, name, binding_method, data);
490 #endif
491 #endif
492
493 slave->trace_id = ((MESH_MEM_INFO *)master->mem_info)->next_trace_id++;
494
495 /****************************************************************************/
496 /* Turn the submesh into a parametric mesh if necessary. */
497 /****************************************************************************/
498 if (master->parametric) {
499 master->parametric->inherit_parametric(slave);
500 }
501
502 return slave;
503 }
504
505 /****************************************************************************/
506 /* get_bndry_submesh(master, name): */
507 /* Convert the entire boundary into a submesh. */
508 /****************************************************************************/
bndry_binding_method(MESH * master,MACRO_EL * mel,int wall,void * data)509 static bool bndry_binding_method(MESH *master,
510 MACRO_EL *mel, int wall, void *data)
511 {
512 if (!mel->neigh[wall]) {
513 return true;
514 } else {
515 return false;
516 }
517 }
518
get_bndry_submesh(MESH * master,const char * name)519 MESH *get_bndry_submesh(MESH *master, const char *name)
520 {
521 return get_submesh(master, name, bndry_binding_method, NULL);
522 }
523
read_bndry_submesh(MESH * master,const char * slave_filename)524 MESH *read_bndry_submesh(MESH *master, const char *slave_filename)
525 {
526 return read_submesh(master, slave_filename, bndry_binding_method, NULL);
527 }
528
read_bndry_submesh_xdr(MESH * master,const char * slave_filename)529 MESH *read_bndry_submesh_xdr(MESH *master, const char *slave_filename)
530 {
531 return read_submesh_xdr(master, slave_filename, bndry_binding_method, NULL);
532 }
533
534 /****************************************************************************/
535 /* get_bndry_submesh_by_type(master, name, type): */
536 /* Convert boundary facets of the given type into a submesh. */
537 /****************************************************************************/
bndry_binding_method_by_type(MESH * master,MACRO_EL * mel,int wall,void * data)538 static bool bndry_binding_method_by_type(MESH *master,
539 MACRO_EL *mel, int wall, void *data)
540 {
541 BNDRY_TYPE type = *(BNDRY_TYPE *)data;
542
543 if (mel->wall_bound[wall] == type) {
544 return true;
545 } else {
546 return false;
547 }
548 }
549
get_bndry_submesh_by_type(MESH * master,const char * name,BNDRY_TYPE type)550 MESH *get_bndry_submesh_by_type(MESH *master, const char *name, BNDRY_TYPE type)
551 {
552 return get_submesh(master,
553 name, bndry_binding_method_by_type, (void *)&type);
554 }
555
read_bndry_submesh_by_type(MESH * master,const char * slave_filename,BNDRY_TYPE type)556 MESH *read_bndry_submesh_by_type(MESH *master,
557 const char *slave_filename, BNDRY_TYPE type)
558 {
559 union
560 {
561 BNDRY_TYPE type;
562 void *data;
563 } voidptype;
564
565 voidptype.type = type;
566
567 return read_submesh(master, slave_filename,
568 bndry_binding_method_by_type,
569 voidptype.data);
570 }
571
read_bndry_submesh_by_type_xdr(MESH * master,const char * slave_filename,BNDRY_TYPE type)572 MESH *read_bndry_submesh_by_type_xdr(MESH *master,
573 const char *slave_filename,
574 BNDRY_TYPE type)
575 {
576 union
577 {
578 BNDRY_TYPE type;
579 void *data;
580 } voidptype;
581
582 voidptype.type = type;
583
584 return read_submesh_xdr(master, slave_filename,
585 bndry_binding_method_by_type,
586 voidptype.data);
587 }
588
589 /******************************************************************************
590 * get_bndry_submesh_by_type(master, name, segment):
591 *
592 * Convert the boundary segments mentioned in the bit-mask segment
593 * into a trace-mesh.
594 ******************************************************************************/
bndry_binding_method_by_segment(MESH * master,MACRO_EL * mel,int wall,void * data)595 static bool bndry_binding_method_by_segment(MESH *master,
596 MACRO_EL *mel, int wall, void *data)
597 {
598 const FLAGS *segment = (const FLAGS *)data;
599
600 if (BNDRY_FLAGS_IS_AT_BNDRY(segment, mel->wall_bound[wall])) {
601 return true;
602 } else {
603 return false;
604 }
605 }
606
get_bndry_submesh_by_segment(MESH * master,const char * name,const BNDRY_FLAGS segment)607 MESH *get_bndry_submesh_by_segment(MESH *master,
608 const char *name, const BNDRY_FLAGS segment)
609 {
610 return get_submesh(master, name,
611 bndry_binding_method_by_segment, (void *)segment);
612 }
613
read_bndry_submesh_by_segment(MESH * master,const char * slave_filename,const BNDRY_FLAGS segment)614 MESH *read_bndry_submesh_by_segment(MESH *master,
615 const char *slave_filename,
616 const BNDRY_FLAGS segment)
617 {
618 return read_submesh(master, slave_filename,
619 bndry_binding_method_by_segment, (void *)segment);
620 }
621
read_bndry_submesh_by_segment_xdr(MESH * master,const char * slave_filename,const BNDRY_FLAGS segment)622 MESH *read_bndry_submesh_by_segment_xdr(MESH *master,
623 const char *slave_filename,
624 const BNDRY_FLAGS segment)
625 {
626 return read_submesh_xdr(master, slave_filename,
627 bndry_binding_method_by_segment, (void *)segment);
628 }
629
630 /****************************************************************************/
631 /* unchain_submesh(slave): This routine destroys the ties between a slave */
632 /* mesh and its master mesh. It DOES NOT destroy the slave mesh! That has to*/
633 /* be done by another call to free_mesh(). After this call, refinement and */
634 /* coarsening of master and slave are independent, and there is no further */
635 /* tie between the two meshes. */
636 /****************************************************************************/
637
unchain_submesh(MESH * slave)638 void unchain_submesh(MESH *slave)
639 {
640 FUNCNAME("unchain_submesh");
641 MESH *master;
642 MESH_MEM_INFO *master_info, *slave_info;
643 int i;
644
645 /****************************************************************************/
646 /* Do the checks for obvious user errors. */
647 /****************************************************************************/
648
649 if (!slave) {
650 ERROR("No slave mesh specified!\n");
651 return;
652 }
653
654 slave_info = (MESH_MEM_INFO *)slave->mem_info;
655
656 if (!(master = slave_info->master)) {
657 ERROR("This mesh is not a slave mesh!\n");
658 return;
659 }
660
661 /****************************************************************************/
662 /* Look for the slave and the master mesh->mem_info. */
663 /****************************************************************************/
664
665 master_info = (MESH_MEM_INFO *)master->mem_info;
666
667 for (i = 0; i < master_info->n_slaves; i++)
668 if (master_info->slaves[i] == slave)
669 break;
670 TEST_EXIT(i < master_info->n_slaves,
671 "Could not find the slave mesh in slave vector!\n");
672
673 /****************************************************************************/
674 /* Call the parametric unchain-hook if necessary */
675 /****************************************************************************/
676 if (slave->parametric && slave->parametric->unchain_parametric)
677 slave->parametric->unchain_parametric(slave);
678
679 /****************************************************************************/
680 /* Resize the slaves vector on the master mesh->mem_info */
681 /****************************************************************************/
682
683 for (; i < master_info->n_slaves - 1; i++)
684 master_info->slaves[i] = master_info->slaves[i+1];
685
686 if (master_info->n_slaves > 1)
687 master_info->slaves = MEM_REALLOC(master_info->slaves,
688 master_info->n_slaves,
689 master_info->n_slaves - 1,
690 MESH *);
691 else {
692 MEM_FREE(master_info->slaves, 1, MESH *);
693
694 master_info->slaves = NULL;
695 }
696 master_info->n_slaves--;
697
698 /****************************************************************************/
699 /* Free the DOF_PTR_VECs and fe_spaces */
700 /****************************************************************************/
701
702 free_dof_ptr_vec(slave_info->master_binding);
703 free_dof_ptr_vec(slave_info->slave_binding);
704
705 /****************************************************************************/
706 /* Set the necessary pointers to NULL, so that this mesh is now free. */
707 /****************************************************************************/
708
709 slave_info->master = NULL;
710 slave_info->master_binding = NULL;
711 slave_info->slave_binding = NULL;
712
713 slave->trace_id = -1;
714
715 return;
716 }
717
718
read_submesh(MESH * master,const char * slave_filename,bool (* binding_method)(MESH * master,MACRO_EL * el,int wall,void * data),void * data)719 MESH *read_submesh(MESH *master,
720 const char *slave_filename,
721 bool (*binding_method)(MESH *master,
722 MACRO_EL *el, int wall,
723 void *data),
724 void *data)
725 {
726 return read_submesh_gen(false, master, slave_filename, binding_method, data);
727 }
728
read_submesh_xdr(MESH * master,const char * slave_filename,bool (* binding_method)(MESH * master,MACRO_EL * el,int wall,void * data),void * data)729 MESH *read_submesh_xdr(MESH *master,
730 const char *slave_filename,
731 bool (*binding_method)(MESH *master,
732 MACRO_EL *el, int wall,
733 void *data),
734 void *data)
735 {
736 return read_submesh_gen(true, master, slave_filename, binding_method, data);
737 }
738
739 /* Return a pointer to an existing trace-mesh with ID if this is
740 * available, NULL otherwise.
741 */
lookup_submesh_by_id(MESH * mesh,int id)742 MESH *lookup_submesh_by_id(MESH *mesh, int id)
743 {
744 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
745 int i;
746
747 for (i = 0; i < mem_info->n_slaves; i++) {
748 if (mem_info->slaves[i]->trace_id == id) {
749 return mem_info->slaves[i];
750 }
751 }
752 return NULL;
753 }
754
lookup_submesh_by_name(MESH * mesh,const char * name)755 MESH *lookup_submesh_by_name(MESH *mesh, const char *name)
756 {
757 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
758 int i;
759
760 for (i = 0; i < mem_info->n_slaves; i++) {
761 if (mem_info->slaves[i]->name == NULL) {
762 continue;
763 }
764 if (strcmp(mem_info->slaves[i]->name, name) == 0) {
765 return mem_info->slaves[i];
766 }
767 }
768 return NULL;
769 }
770
771 /* Search for an existing trace mesh conforming to BINDING_METHOD(). */
lookup_submesh_by_binding(MESH * master,bool (* binding_method)(MESH * master,MACRO_EL * el,int wall,void * data),void * data)772 MESH *lookup_submesh_by_binding(MESH *master,
773 bool (*binding_method)(MESH *master,
774 MACRO_EL *el, int wall,
775 void *data),
776 void *data)
777 {
778 MESH_MEM_INFO *m_info = (MESH_MEM_INFO *)master->mem_info;
779 int m_dim = master->dim;
780 int k;
781
782 for (k = 0; k < m_info->n_slaves; k++) {
783 MESH *slave = m_info->slaves[k];
784 MACRO_EL *s_mel = slave->macro_els;
785 MACRO_EL *s_mel_end = s_mel + slave->n_macro_el;
786 int n;
787 bool ok = true;
788
789 for (n = 0; ok && n < master->n_macro_el; ++n) {
790 MACRO_EL *m_mel = master->macro_els + n;
791 int i;
792 for (i = 0; i < N_NEIGH(m_dim); i++) {
793 if (binding_method(master, m_mel, i, data)) {
794 if (s_mel >= s_mel_end ||
795 s_mel->master.macro_el != m_mel ||
796 s_mel->master.opp_vertex != i) {
797 ok = false;
798 }
799 s_mel++;
800 }
801 }
802 }
803 if (ok == true && s_mel == s_mel_end) {
804 /* At this point all slave macro element have found their
805 * matching master element, and no (master, wall) subject to
806 * binding_method() is without slave.
807 */
808 return slave;
809 }
810 }
811 return NULL;
812 }
813
lookup_bndry_submesh_by_type(MESH * master,BNDRY_TYPE type)814 MESH *lookup_bndry_submesh_by_type(MESH *master, BNDRY_TYPE type)
815 {
816 return lookup_submesh_by_binding(master,
817 bndry_binding_method_by_type, (void *)&type);
818 }
819
lookup_bndry_submesh_by_segment(MESH * master,const BNDRY_FLAGS segment)820 MESH *lookup_bndry_submesh_by_segment(MESH *master, const BNDRY_FLAGS segment)
821 {
822 return lookup_submesh_by_binding(master,
823 bndry_binding_method_by_segment,
824 (void *)segment);
825 }
826
lookup_bndry_submesh(MESH * master)827 MESH *lookup_bndry_submesh(MESH *master)
828 {
829 return lookup_submesh_by_binding(master, bndry_binding_method, NULL);
830 }
831
832 #if 0
833 /****************************************************************************/
834 /* trace_dof_real_vec[_d](slave_vec, master_vec): */
835 /* A routine for transporting the values of "master_vec" at DOFs along the */
836 /* interface to a slave mesh to "slave_vec". This could be described as a */
837 /* discrete trace operation. */
838 /****************************************************************************/
839
840 void trace_dof_real_vec_old(DOF_REAL_VEC *slave_vec, DOF_REAL_VEC *master_vec)
841 {
842 FUNCNAME("trace_dof_real_vec");
843
844 MESH *slave = NULL, *master = NULL;
845 TRAVERSE_STACK *m_stack = get_traverse_stack();
846 EL_INFO s_el_info = { NULL, };
847 DOF_PTR_VEC *s_dpv, *m_dpv;
848 const DOF_ADMIN *s_admin, *s_ptr_admin, *m_ptr_admin;
849 const BAS_FCTS *s_bfcts, *m_bfcts;
850 int s_n_bas_fcts, m_dim;
851
852 /****************************************************************************/
853 /* Check for user errors. */
854 /****************************************************************************/
855 TEST_EXIT(slave_vec,"No slave DOF_REAL_VEC given!\n");
856 TEST_EXIT(m_vec = master_vec,"No master DOF_REAL_VEC given!\n");
857
858 /****************************************************************************/
859 /* Set the necessary pointers and allocate two TRAVERSE_STACKS. */
860 /****************************************************************************/
861 slave = slave_vec->fe_space->mesh;
862 master = master_vec->fe_space->mesh;
863
864 TEST_EXIT(((MESH_MEM_INFO *)slave->mem_info)->master == master,
865 "Master and slave vectors do not seem to fit together!\n");
866
867 m_dim = master->dim;
868
869 s_dpv = ((MESH_MEM_INFO *)slave->mem_info)->master_binding;
870 m_dpv = ((MESH_MEM_INFO *)slave->mem_info)->slave_binding;
871
872 s_admin = slave_vec->fe_space->admin;
873 s_bfcts = slave_vec->fe_space->bas_fcts;
874 m_bfcts = master_vec->fe_space->bas_fcts;
875 s_ptr_admin = s_dpv->fe_space->admin;
876 m_ptr_admin = m_dpv->fe_space->admin;
877
878 s_n_bas_fcts = s_bfcts->n_bas_fcts;
879
880 TEST_EXIT(s_bfcts->interpol,
881 "Interpolation routine not present for slave vector!\n");
882
883 TEST_EXIT(s_bfcts->get_dof_indices,
884 "'get_dof_indices' routine not present for slave vector!\n");
885
886 TEST_EXIT(get_local_m_vec = master_vec->fe_space->bas_fcts->get_real_vec,
887 "'get_real_vec' routine not present for master vector!\n");
888
889 m_el_info = traverse_first(m_stack, master, -1,
890 CALL_LEAF_EL|FILL_ORIENTATION);
891 s_el_info.mesh = slave;
892
893 {
894 int i;
895 REAL s_local_vec[s_n_bas_fcts];
896 DOF s_local_dofs[s_n_bas_fcts];
897
898 while(m_el_info) {
899
900 if (INIT_ELEMENT(m_el_info, m_bfcts) == INIT_EL_TAG_NULL) {
901 continue;
902 }
903
904 for (m_subsimplex = 0; m_subsimplex < N_NEIGH(m_dim); m_subsimplex++) {
905 /****************************************************************************/
906 /* Look for a slave element on edge/face m_subsimplex. */
907 /****************************************************************************/
908 switch(m_dim) {
909 case 1:
910 if ((s_el_info.el = (EL *)
911 m_dpv->vec[m_el_info->el->dof[master->node[VERTEX] + m_subsimplex]
912 [m_ptr_admin->n0_dof[VERTEX]]]))
913 break;
914 case 2:
915 if ((s_el_info.el = (EL *)
916 m_dpv->vec[m_el_info->el->dof[master->node[EDGE] + m_subsimplex]
917 [m_ptr_admin->n0_dof[EDGE]]]))
918 goto found_slave_element;
919 break;
920 case 3:
921 if ((s_el_info.el = (EL *)
922 m_dpv->vec[m_el_info->el->dof[master->node[FACE] + m_subsimplex]
923 [m_ptr_admin->n0_dof[FACE]]]))
924 goto found_slave_element;
925 break;
926 }
927 continue;
928
929 found_slave_element:
930 /****************************************************************************/
931 /* Interpolate the master vector into s_local_vec and copy these values */
932 /* to the correct position in slave_vec->vec. */
933 /****************************************************************************/
934 GET_DOF_INDICES(s_bfcts, s_el_info.el, s_admin, s_local_dofs);
935 INTERPOL(s_bfcts,
936 &s_el_info, 0, NULL, NULL, trace_loc, NULL, s_local_vec);
937
938 for (i = 0; i < s_n_bas_fcts; i++)
939 slave_vec->vec[s_local_dofs[i]] = s_local_vec[i];
940 }
941
942 /****************************************************************************/
943 /* Advance the slave stack. */
944 /****************************************************************************/
945 m_el_info = traverse_next(m_stack, m_el_info);
946 }
947 }
948
949
950 /****************************************************************************/
951 /* Clean up. */
952 /****************************************************************************/
953 free_traverse_stack(m_stack);
954
955 return;
956 }
957
958
959 void trace_dof_real_d_vec_old(DOF_REAL_D_VEC *slave_vec,
960 DOF_REAL_D_VEC *master_vec)
961 {
962 FUNCNAME("trace_dof_real_d_vec");
963
964 MESH *slave = NULL, *master = NULL;
965 TRAVERSE_STACK *m_stack = get_traverse_stack();
966 EL_INFO s_el_info = { NULL, };
967 DOF_PTR_VEC *s_dpv, *m_dpv;
968 const DOF_ADMIN *s_admin, *s_ptr_admin, *m_ptr_admin;
969 const BAS_FCTS *s_bfcts, *m_bfcts;
970 int s_n_bas_fcts, m_dim;
971
972 /****************************************************************************/
973 /* Check for user errors. */
974 /****************************************************************************/
975 TEST_EXIT(slave_vec,"No slave DOF_REAL_D_VEC given!\n");
976 TEST_EXIT(m_vec_d = master_vec,"No master DOF_REAL_D_VEC given!\n");
977
978 /****************************************************************************/
979 /* Set the necessary pointers and allocate two TRAVERSE_STACKS. */
980 /****************************************************************************/
981 slave = slave_vec->fe_space->mesh;
982 master = master_vec->fe_space->mesh;
983
984 TEST_EXIT(((MESH_MEM_INFO *)slave->mem_info)->master == master,
985 "Master and slave vectors do not seem to fit together!\n");
986
987 m_dim = master->dim;
988
989 s_dpv = ((MESH_MEM_INFO *)slave->mem_info)->master_binding;
990 m_dpv = ((MESH_MEM_INFO *)slave->mem_info)->slave_binding;
991
992 s_admin = slave_vec->fe_space->admin;
993 s_bfcts = slave_vec->fe_space->bas_fcts;
994 s_bfcts = master_vec->fe_space->bas_fcts;
995 m_bfcts = master_vec->fe_space->bas_fcts;
996 s_ptr_admin = s_dpv->fe_space->admin;
997 m_ptr_admin = m_dpv->fe_space->admin;
998
999 s_n_bas_fcts = slave_vec->fe_space->bas_fcts->n_bas_fcts;
1000
1001 TEST_EXIT(s_bfcts->interpol_d,
1002 "Interpolation routine not present for slave vector!\n");
1003
1004 TEST_EXIT(s_bfcts->get_dof_indices,
1005 "'get_dof_indices' routine not present for slave vector!\n");
1006
1007 TEST_EXIT(get_local_m_vec_d = master_vec->fe_space->bas_fcts->get_real_d_vec,
1008 "'get_real_d_vec' routine not present for master vector!\n");
1009
1010 m_el_info = traverse_first(m_stack, master, -1,
1011 CALL_LEAF_EL|FILL_ORIENTATION);
1012 s_el_info.mesh = slave;
1013
1014 {
1015 int i;
1016 REAL_D s_local_vec_d[s_n_bas_fcts];
1017 DOF s_local_dofs[s_n_bas_fcts];
1018
1019 while(m_el_info) {
1020
1021 if (INIT_ELEMENT(m_el_info, m_bfcts) == INIT_EL_TAG_NULL) {
1022 continue;
1023 }
1024
1025 for (m_subsimplex = 0; m_subsimplex < N_NEIGH(m_dim); m_subsimplex++) {
1026 /****************************************************************************/
1027 /* Look for a slave element on edge/face m_subsimplex. */
1028 /****************************************************************************/
1029 switch(m_dim) {
1030 case 1:
1031 if ((s_el_info.el = (EL *)
1032 m_dpv->vec[m_el_info->el->dof[master->node[VERTEX] + m_subsimplex]
1033 [m_ptr_admin->n0_dof[VERTEX]]]))
1034 break;
1035 case 2:
1036 if ((s_el_info.el = (EL *)
1037 m_dpv->vec[m_el_info->el->dof[master->node[EDGE] + m_subsimplex]
1038 [m_ptr_admin->n0_dof[EDGE]]]))
1039 goto found_slave_element;
1040 break;
1041 case 3:
1042 if ((s_el_info.el = (EL *)
1043 m_dpv->vec[m_el_info->el->dof[master->node[FACE] + m_subsimplex]
1044 [m_ptr_admin->n0_dof[FACE]]]))
1045 goto found_slave_element;
1046 break;
1047 }
1048 continue;
1049
1050 found_slave_element:
1051 /****************************************************************************/
1052 /* Interpolate the master vector into s_local_vec and copy these values */
1053 /* to the correct position in slave_vec->vec. */
1054 /****************************************************************************/
1055 GET_DOF_INDICES(s_bfcts, s_el_info.el, s_admin, s_local_dofs);
1056 INTERPOL_D(s_bfcts,
1057 &s_el_info, 0, NULL, NULL, trace_loc_d, NULL, s_local_vec_d);
1058
1059 for (i = 0; i < s_n_bas_fcts; i++)
1060 COPY_DOW(s_local_vec_d[i], slave_vec->vec[s_local_dofs[i]]);
1061 }
1062
1063 /****************************************************************************/
1064 /* Advance the slave stack. */
1065 /****************************************************************************/
1066 m_el_info = traverse_next(m_stack, m_el_info);
1067 }
1068 }
1069
1070 /****************************************************************************/
1071 /* Clean up. */
1072 /****************************************************************************/
1073 free_traverse_stack(m_stack);
1074
1075 return;
1076 }
1077
1078 #endif
1079
1080 /****************************************************************************/
1081 /* get_slave_dof_mapping(m_fe_space, s_map): Fill s_map on the slave mesh */
1082 /* with the corresponding DOFs of m_fe_space. This is only implemented for */
1083 /* Lagrange FE spaces. Master and slave fe_spaces must have the same degree.*/
1084 /* */
1085 /* s_map is not updated during mesh or DOF (dof_compress(master)) changes!! */
1086 /****************************************************************************/
1087
get_slave_dof_mapping(const FE_SPACE * m_fe_space,DOF_INT_VEC * s_map)1088 void get_slave_dof_mapping(const FE_SPACE *m_fe_space, DOF_INT_VEC *s_map)
1089 {
1090 FUNCNAME("get_slave_dof_mapping");
1091 MESH *master;
1092 MESH *slave;
1093 const FE_SPACE *s_fe_space;
1094 const BAS_FCTS *m_bfcts, *s_bfcts;
1095 const DOF_ADMIN *m_admin, *s_admin;
1096 DOF_PTR_VEC *m_dpv, *s_dpv;
1097 TRAVERSE_STACK *m_stack = get_traverse_stack();
1098 const EL_INFO *m_el_info;
1099 DOF *m_local_dofs, *s_local_dofs;
1100 int m_dim, m_n_dofs, s_n_dofs, m_subsimplex;
1101 int m_n0, s_n0, m_n, s_n;
1102 const EL *s_el, *m_el;
1103 FLAGS which_elements;
1104
1105 /****************************************************************************/
1106 /* Check for user errors. */
1107 /****************************************************************************/
1108 TEST_EXIT(m_fe_space,"No master FE_SPACE given!\n");
1109 TEST_EXIT(s_map,"No DOF_INT_VEC s_map given!\n");
1110
1111 s_fe_space = s_map->fe_space;
1112
1113 TEST_EXIT(s_fe_space,"No slave FE_SPACE found!\n");
1114
1115 m_admin = m_fe_space->admin;
1116 s_admin = s_fe_space->admin;
1117
1118 m_bfcts = m_fe_space->bas_fcts;
1119 s_bfcts = s_fe_space->bas_fcts;
1120
1121 TEST_EXIT(m_bfcts,
1122 "Sorry, only implemented for FE_SPACEs derived from local "
1123 "basis functions\n");
1124
1125 TEST_EXIT(s_bfcts == m_bfcts->trace_bas_fcts,
1126 "s_map->fe_space->bas_fcts != m_bfcts->trace_bas_fcts.\n");
1127
1128 master = m_fe_space->mesh;
1129 slave = s_fe_space->mesh;
1130 m_dim = master->dim;
1131
1132 TEST_EXIT(((MESH_MEM_INFO *)slave->mem_info)->master == master,
1133 "Master and slave meshes do not seem to belong together!\n");
1134 TEST_EXIT(strstr(m_fe_space->bas_fcts->name, "lagrange")
1135 && strstr(s_fe_space->bas_fcts->name, "lagrange"),
1136 "Sorry, only implemented for Lagrange Finite Elements!\n");
1137
1138 TEST_EXIT(m_admin->flags == s_admin->flags, "different flag values!\n");
1139 if (s_admin->flags & ADM_PRESERVE_COARSE_DOFS)
1140 which_elements = CALL_EVERY_EL_PREORDER;
1141 else
1142 which_elements = CALL_LEAF_EL;
1143
1144 /****************************************************************************/
1145 /* Initialize values. */
1146 /****************************************************************************/
1147
1148 /* Mark everything as unused. */
1149 FOR_ALL_DOFS(s_admin, s_map->vec[dof] = -1);
1150
1151 m_dpv = ((MESH_MEM_INFO *)slave->mem_info)->slave_binding;
1152 s_dpv = ((MESH_MEM_INFO *)slave->mem_info)->master_binding;
1153
1154 s_n0 = s_dpv->fe_space->admin->n0_dof[CENTER];
1155 s_n = slave->node[CENTER];
1156
1157 m_n_dofs = m_fe_space->bas_fcts->n_bas_fcts;
1158 s_n_dofs = s_fe_space->bas_fcts->n_bas_fcts;
1159
1160 m_local_dofs = MEM_ALLOC(m_n_dofs, DOF);
1161 s_local_dofs = MEM_ALLOC(s_n_dofs, DOF);
1162
1163 switch(m_dim) {
1164 case 1:
1165 m_n0 = m_dpv->fe_space->admin->n0_dof[VERTEX];
1166 m_n = master->node[VERTEX];
1167
1168 for (m_el_info = traverse_first(m_stack, master, -1, which_elements);
1169 m_el_info;
1170 m_el_info = traverse_next(m_stack, m_el_info)) {
1171
1172 if (INIT_ELEMENT(m_el_info, m_bfcts) == INIT_EL_TAG_NULL) {
1173 continue;
1174 }
1175
1176 m_el = m_el_info->el;
1177 GET_DOF_INDICES(m_bfcts, m_el, m_admin, m_local_dofs);
1178
1179 for (m_subsimplex = 0; m_subsimplex < N_NEIGH_1D; m_subsimplex++) {
1180 s_el = (EL *)m_dpv->vec[m_el->dof[m_n + m_subsimplex][m_n0]];
1181 if (s_el && m_el==(EL *)s_dpv->vec[s_el->dof[s_n][s_n0]]) {
1182 GET_DOF_INDICES(s_bfcts, s_el, s_admin, s_local_dofs);
1183 s_map->vec[s_local_dofs[0]] = m_local_dofs[m_subsimplex];
1184 }
1185 }
1186 }
1187 break;
1188 #if DIM_MAX > 1
1189 case 2:
1190 m_n0 = m_dpv->fe_space->admin->n0_dof[EDGE];
1191 m_n = master->node[EDGE];
1192
1193 for (m_el_info = traverse_first(m_stack, master, -1, which_elements);
1194 m_el_info;
1195 m_el_info = traverse_next(m_stack, m_el_info)) {
1196
1197 if (INIT_ELEMENT(m_el_info, m_bfcts) == INIT_EL_TAG_NULL) {
1198 continue;
1199 }
1200
1201 m_el = m_el_info->el;
1202 GET_DOF_INDICES(m_bfcts, m_el, m_admin, m_local_dofs);
1203
1204 for (m_subsimplex = 0; m_subsimplex < N_NEIGH_2D; m_subsimplex++) {
1205 int i;
1206
1207 s_el = (EL *)m_dpv->vec[m_el->dof[m_n + m_subsimplex][m_n0]];
1208 if (s_el && m_el==(EL *)s_dpv->vec[s_el->dof[s_n][s_n0]]) {
1209 GET_DOF_INDICES(s_bfcts, s_el, s_admin, s_local_dofs);
1210
1211 for (i = 0; i < s_n_dofs; i++)
1212 s_map->vec[s_local_dofs[i]] =
1213 m_local_dofs[m_bfcts->trace_dof_map[0][0][m_subsimplex][i]];
1214 }
1215 }
1216 }
1217 break;
1218 #if DIM_MAX > 2
1219 case 3:
1220 m_n0 = m_dpv->fe_space->admin->n0_dof[FACE];
1221 m_n = master->node[FACE];
1222
1223 for (m_el_info = traverse_first(m_stack, master, -1,
1224 which_elements|FILL_ORIENTATION);
1225 m_el_info;
1226 m_el_info = traverse_next(m_stack, m_el_info)) {
1227 int type, i;
1228
1229 if (INIT_ELEMENT(m_el_info, m_bfcts) == INIT_EL_TAG_NULL) {
1230 continue;
1231 }
1232
1233 m_el = m_el_info->el;
1234 GET_DOF_INDICES(m_bfcts, m_el, m_admin, m_local_dofs);
1235
1236 if (m_el_info->el_type == 0)
1237 type = 0;
1238 else
1239 type = 1;
1240
1241 for (m_subsimplex = 0; m_subsimplex < N_NEIGH_3D; m_subsimplex++) {
1242 s_el = (EL *)m_dpv->vec[m_el->dof[m_n + m_subsimplex][m_n0]];
1243 if (s_el && m_el==(EL *)s_dpv->vec[s_el->dof[s_n][s_n0]]) {
1244 GET_DOF_INDICES(s_bfcts, s_el, s_admin, s_local_dofs);
1245
1246 for (i = 0; i < s_n_dofs; i++) {
1247 if (m_el_info->orientation > 0)
1248 s_map->vec[s_local_dofs[i]] =
1249 m_local_dofs[
1250 m_bfcts->trace_dof_map[type][0][m_subsimplex][i]];
1251 else
1252 s_map->vec[s_local_dofs[i]] =
1253 m_local_dofs[
1254 m_bfcts->trace_dof_map[type][1][m_subsimplex][i]];
1255 }
1256 }
1257 }
1258 }
1259 break;
1260 #endif
1261 #endif
1262 default:
1263 ERROR_EXIT("Illegal dimension!\n");
1264 }
1265
1266 /****************************************************************************/
1267 /* Clean up and return. */
1268 /****************************************************************************/
1269 free_traverse_stack(m_stack);
1270 MEM_FREE(m_local_dofs, m_n_dofs, DOF);
1271 MEM_FREE(s_local_dofs, s_n_dofs, DOF);
1272
1273 return;
1274 }
1275
1276 /****************************************************************************/
1277 /* get_master(slave): Return the master of slave if present. */
1278 /****************************************************************************/
1279
get_master(MESH * slave)1280 MESH *get_master(MESH *slave)
1281 {
1282 return ((MESH_MEM_INFO *)slave->mem_info)->master;
1283 }
1284
get_slave_el(const EL * el,int wall,MESH * trace_mesh)1285 const EL *get_slave_el(const EL *el, int wall, MESH *trace_mesh)
1286 {
1287 MESH_MEM_INFO *trace_minfo = trace_mesh->mem_info;
1288 MESH *mesh = get_master(trace_mesh);
1289 const DOF_PTR_VEC *trace_binding = trace_minfo->slave_binding;
1290 int n0, node;
1291
1292 switch (trace_mesh->dim) {
1293 case 2: node = FACE; break;
1294 case 1: node = EDGE; break;
1295 case 0: node = VERTEX; break;
1296 default: node = -1; break;
1297 }
1298
1299 trace_binding = ((MESH_MEM_INFO *)trace_mesh->mem_info)->slave_binding;
1300 n0 = trace_binding->fe_space->admin->n0_dof[node];
1301 node = mesh->node[node];
1302
1303 return (const EL *)trace_binding->vec[el->dof[node+wall][n0]];
1304 }
1305
fill_slave_el_info(EL_INFO * slv_el_info,const EL_INFO * el_info,int wall,MESH * trace_mesh)1306 void fill_slave_el_info(EL_INFO *slv_el_info,
1307 const EL_INFO *el_info, int wall, MESH *trace_mesh)
1308 {
1309 const EL *slv_el = get_slave_el(el_info->el, wall, trace_mesh);
1310 int dim = trace_mesh->dim;
1311
1312 slv_el_info->fill_flag = FILL_NOTHING;
1313 slv_el_info->mesh = trace_mesh;
1314
1315 slv_el_info->macro_el = NULL; /* FIXME */
1316 slv_el_info->el = (EL *)slv_el;
1317 slv_el_info->parent = NULL /* could do better here */;
1318
1319 slv_el_info->master.el = el_info->el;
1320 slv_el_info->master.opp_vertex = wall;
1321 slv_el_info->master.el_type = el_info->el_type;
1322 slv_el_info->master.orientation = el_info->orientation;
1323 slv_el_info->fill_flag |= FILL_MASTER_INFO;
1324
1325 slv_el_info->el_geom_cache.fill_flag = 0U;
1326 slv_el_info->el_geom_cache.current_el = slv_el_info->el;
1327
1328 if (el_info->fill_flag & FILL_COORDS) {
1329 if (dim == 2) {
1330 #if DIM_MAX > 2
1331 int mst_t, mst_o, v;
1332
1333 mst_t = el_info->el_type;
1334 mst_o = el_info->orientation;
1335
1336 for (v = 0; v < N_VERTICES_3D; v++) {
1337 int sv = slave_numbering_3d[!!mst_t][mst_o < 0][wall][v];
1338 if (sv >= 0) {
1339 COPY_DOW(el_info->coord[v], slv_el_info->coord[sv]);
1340 }
1341 }
1342 #endif
1343 } else {
1344 int v;
1345 for (v = 0; v < N_VERTICES(dim); v++) {
1346 int mv = (wall + v) % N_VERTICES(dim+1);
1347 COPY_DOW(el_info->coord[mv], slv_el_info->coord[v]);
1348 }
1349 }
1350 COPY_DOW(el_info->coord[wall], slv_el_info->master.opp_coord);
1351 slv_el_info->fill_flag |= FILL_COORDS;
1352 }
1353 if (el_info->fill_flag & FILL_NEIGH) {
1354 slv_el_info->mst_neigh.el = el_info->neigh[wall];
1355 if (el_info->neigh[wall]) {
1356 slv_el_info->mst_neigh.opp_vertex = el_info->opp_vertex[wall];
1357 slv_el_info->mst_neigh.el_type = 0;
1358 slv_el_info->mst_neigh.orientation = 1;
1359 if (el_info->fill_flag & FILL_OPP_COORDS) {
1360 COPY_DOW(el_info->opp_coord[wall], slv_el_info->mst_neigh.opp_coord);
1361 }
1362 }
1363 slv_el_info->fill_flag |= FILL_MASTER_NEIGH;
1364 }
1365 }
1366
1367 /* el_info has to refer to a trace element, lambda are trace
1368 * coordinates.
1369 */
trace_to_bulk_coords_2d(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)1370 void trace_to_bulk_coords_2d(REAL_B result,
1371 const REAL_B lambda,
1372 const EL_INFO *el_info)
1373 {
1374 #if DIM_MAX > 2
1375 int mst_ov = el_info->master.opp_vertex;
1376 int mst_t = el_info->master.el_type;
1377 int mst_o = el_info->master.orientation;
1378 int v;
1379
1380 for (v = 0; v < N_LAMBDA_3D; v++) {
1381 int sv = slave_numbering_3d[!!mst_t][mst_o < 0][mst_ov][v];
1382 if (sv >= 0) {
1383 result[v] = lambda[sv];
1384 }
1385 }
1386 result[mst_ov] = 0.0;
1387
1388 #else
1389 FUNCNAME("trace_to_bulk_coords_2d");
1390 ERROR_EXIT("This must not happen.\n");
1391 #endif
1392 }
1393
1394 /* el_info has to refer to a trace element, lambda are bulk
1395 * coordinates.
1396 */
bulk_to_trace_coords_2d(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)1397 void bulk_to_trace_coords_2d(REAL_B result,
1398 const REAL_B lambda,
1399 const EL_INFO *el_info)
1400 {
1401 FUNCNAME("bulk_to_bulk_coords_2d");
1402 #if DIM_MAX > 2
1403 int mst_ov = el_info->master.opp_vertex;
1404 int mst_t = el_info->master.el_type;
1405 int mst_o = el_info->master.orientation;
1406 int v;
1407
1408 DEBUG_TEST_EXIT(lambda[mst_ov] == 0.0,
1409 "This bulk coordinate does not live on a face.");
1410
1411 for (v = 0; v < N_LAMBDA_3D; v++) {
1412 int sv = slave_numbering_3d[!!mst_t][mst_o < 0][mst_ov][v];
1413 if (sv >= 0) {
1414 result[sv] = lambda[v];
1415 }
1416 }
1417
1418 for (v = N_LAMBDA_2D; v < N_LAMBDA_MAX; v++) {
1419 result[v] = 0.0;
1420 }
1421 #else
1422 ERROR_EXIT("This must not happen.\n");
1423 #endif
1424 }
1425
trace_to_bulk_coords_1d(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)1426 void trace_to_bulk_coords_1d(REAL_B result,
1427 const REAL_B lambda,
1428 const EL_INFO *el_info)
1429 {
1430 #if DIM_MAX > 1
1431 int mst_ov = el_info->master.opp_vertex;
1432 int v;
1433 for (v = 0; v < N_VERTICES_1D; v++) {
1434 int mv = (mst_ov + v + 1) % N_VERTICES_2D;
1435 result[mv] = lambda[v];
1436 }
1437 result[mst_ov] = 0.0;
1438 for (v = N_LAMBDA_2D; v < N_LAMBDA_MAX; v++) {
1439 result[v] = 0.0;
1440 }
1441 #else
1442 FUNCNAME("trace_to_bulk_coords_1d");
1443 ERROR_EXIT("This must not happen.\n");
1444 #endif
1445 }
1446
bulk_to_trace_coords_1d(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)1447 void bulk_to_trace_coords_1d(REAL_B result,
1448 const REAL_B lambda,
1449 const EL_INFO *el_info)
1450 {
1451 FUNCNAME("bulk_to_trace_coords_1d");
1452 #if DIM_MAX > 1
1453 int mst_ov = el_info->master.opp_vertex;
1454 int v;
1455
1456 DEBUG_TEST_EXIT(lambda[mst_ov] == 0.0,
1457 "This bulk coordinate does not live on a face.");
1458
1459 for (v = 0; v < N_VERTICES_1D; v++) {
1460 int mv = (mst_ov + v + 1) % N_VERTICES_2D;
1461 result[v] = lambda[mv];
1462 }
1463
1464 for (v = N_LAMBDA_1D; v < N_LAMBDA_MAX; v++) {
1465 result[v] = 0.0;
1466 }
1467 #else
1468 ERROR_EXIT("This must not happen.\n");
1469 #endif
1470 }
1471
trace_to_bulk_coords_0d(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)1472 void trace_to_bulk_coords_0d(REAL_B result,
1473 const REAL_B lambda,
1474 const EL_INFO *el_info)
1475 {
1476 #if DIM_MAX > 0
1477 int mst_ov = el_info->master.opp_vertex;
1478 int v;
1479
1480 result[mst_ov] = 0.0;
1481 result[1 - mst_ov] = 1.0;
1482 for (v = N_LAMBDA_1D; v < N_LAMBDA_MAX; v++) {
1483 result[v] = 0.0;
1484 }
1485 #else
1486 FUNCNAME("trace_to_bulk_coords_0d");
1487 ERROR_EXIT("This must not happen.\n");
1488 #endif
1489 }
1490
bulk_to_trace_coords_0d(REAL_B result,const REAL_B lambda,const EL_INFO * el_info)1491 void bulk_to_trace_coords_0d(REAL_B result,
1492 const REAL_B lambda,
1493 const EL_INFO *el_info)
1494 {
1495 FUNCNAME("bulk_to_trace_coords_0d");
1496 #if DIM_MAX > 0
1497 int v;
1498
1499 DEBUG_TEST_EXIT(lambda[el_info->master.opp_vertex] == 0.0,
1500 "This bulk coordinate does not live on a face.");
1501
1502 result[0] = 1.0;
1503
1504 for (v = N_LAMBDA_1D; v < N_LAMBDA_MAX; v++) {
1505 result[v] = 0.0;
1506 }
1507 #else
1508 ERROR_EXIT("This must not happen.\n");
1509 #endif
1510 }
1511
1512 /* Construct an EL_INFO structure for the master element attached to
1513 * the given slave element. This requires FILL_MASTER_INFO to be set
1514 * in the fill-flags of el_info.
1515 *
1516 * The master element will have its el_type and orientation set. If
1517 * the slave EL_INFO has co-ordinate information, then also the master
1518 * EL_INFO will be equipped with co-ordinate information.
1519 */
fill_master_el_info(EL_INFO * mst_el_info,const EL_INFO * el_info,FLAGS fill_flags)1520 void fill_master_el_info(EL_INFO *mst_el_info,
1521 const EL_INFO *el_info,
1522 FLAGS fill_flags)
1523 {
1524 FUNCNAME("fill_master_el_info");
1525 int dim = el_info->mesh->dim;
1526 int mst_ov, mst_t = 0, mst_o = 0;
1527
1528 DEBUG_TEST_EXIT(el_info->fill_flag & FILL_MASTER_INFO,
1529 "Master element link not present in "
1530 "slave element descriptor.\n");
1531
1532 memset(mst_el_info, 0, sizeof(*mst_el_info));
1533
1534 mst_el_info->fill_flag = FILL_NOTHING;
1535 mst_el_info->mesh = get_master(el_info->mesh);
1536
1537 mst_el_info->macro_el = el_info->macro_el->master.macro_el;
1538 mst_el_info->el = el_info->master.el;
1539 mst_el_info->parent = NULL /* could do better here */;
1540
1541 mst_el_info->el_geom_cache.fill_flag = 0U;
1542 mst_el_info->el_geom_cache.current_el = mst_el_info->el;
1543
1544 mst_ov = el_info->master.opp_vertex;
1545
1546 #if DIM_MAX > 2
1547 if (dim == 2) {
1548 mst_t = el_info->master.el_type;
1549 mst_o = el_info->master.orientation;
1550
1551 mst_el_info->fill_flag |= FILL_ORIENTATION;
1552
1553 if (fill_flags & FILL_COORDS) {
1554 int v;
1555 for (v = 0; v < N_VERTICES_3D; v++) {
1556 int sv = slave_numbering_3d[!!mst_t][mst_o < 0][mst_ov][v];
1557 if (sv >= 0) {
1558 COPY_DOW(el_info->coord[sv], mst_el_info->coord[v]);
1559 }
1560 }
1561 }
1562 if (fill_flags & FILL_BOUND) {
1563 int v, e;
1564 for (v = 0; v < N_VERTICES_3D; v++) {
1565 int sv = slave_numbering_3d[!!mst_t][mst_o < 0][mst_ov][v];
1566 if (sv >= 0) {
1567 BNDRY_FLAGS_CPY(mst_el_info->vertex_bound[v],
1568 el_info->vertex_bound[sv]);
1569 } else {
1570 BNDRY_FLAGS_INIT(mst_el_info->vertex_bound[v]);
1571 }
1572 }
1573 for (e = 0; e < N_EDGES_3D; e++) {
1574 BNDRY_FLAGS_INIT(mst_el_info->edge_bound[e]);
1575 }
1576 for (e = 0; e < N_EDGES_2D; e++) {
1577 int mst_edge = master_edge_3d[!!mst_t][mst_o < 0][mst_ov][e];
1578 BNDRY_FLAGS_CPY(mst_el_info->edge_bound[mst_edge],
1579 el_info->edge_bound[e]);
1580 }
1581 mst_el_info->face_bound[mst_ov] = el_info->face_bound[0];
1582 mst_el_info->wall_bound[mst_ov] =
1583 el_info->macro_el->master.macro_el->wall_bound[
1584 el_info->macro_el->master.opp_vertex];
1585
1586 mst_el_info->fill_flag |= FILL_BOUND; /* not quite correct, but so what */
1587 }
1588 } else
1589 #endif
1590 {
1591 if (fill_flags & FILL_COORDS) {
1592 int v;
1593 for (v = 0; v < N_VERTICES(dim); v++) {
1594 int mv = (mst_ov + v + 1) % N_VERTICES(dim+1);
1595 COPY_DOW(el_info->coord[v], mst_el_info->coord[mv]);
1596 }
1597 }
1598 if (fill_flags & FILL_BOUND) {
1599 int v;
1600 for (v = 0; v < N_VERTICES(dim); v++) {
1601 int mv = (mst_ov + v + 1) % N_VERTICES(dim+1);
1602 BNDRY_FLAGS_CPY(mst_el_info->vertex_bound[mv],
1603 el_info->vertex_bound[v]);
1604 }
1605 BNDRY_FLAGS_INIT(mst_el_info->vertex_bound[mst_ov]);
1606 if (dim == 1) {
1607 BNDRY_FLAGS_CPY(mst_el_info->edge_bound[mst_ov],
1608 el_info->edge_bound[0]);
1609 }
1610 mst_el_info->wall_bound[mst_ov] =
1611 el_info->macro_el->master.macro_el->wall_bound[
1612 el_info->macro_el->master.opp_vertex];
1613 mst_el_info->fill_flag |= FILL_BOUND; /* not quite correct, but so what */
1614 }
1615 }
1616
1617 if (fill_flags & FILL_NEIGH) {
1618 mst_el_info->neigh[mst_ov] = el_info->mst_neigh.el;
1619 mst_el_info->opp_vertex[mst_ov] = el_info->mst_neigh.opp_vertex;
1620 mst_el_info->fill_flag |= FILL_NEIGH;
1621 if (fill_flags & FILL_OPP_COORDS) {
1622 COPY_DOW(el_info->mst_neigh.opp_coord, mst_el_info->opp_coord[mst_ov]);
1623 mst_el_info->fill_flag |= FILL_OPP_COORDS;
1624 }
1625 }
1626
1627 mst_el_info->el_type = mst_t;
1628 mst_el_info->orientation = mst_o;
1629
1630 if (fill_flags & FILL_COORDS) {
1631 mst_el_info->fill_flag |= FILL_COORDS;
1632 COPY_DOW(el_info->master.opp_coord, mst_el_info->coord[mst_ov]);
1633 }
1634 }
1635
1636 /* Compute the local mapping between global master DOFs and local
1637 * slave DOFs, return value is an array with the global DOFs w.r.t. to
1638 * the master for the local basis functions on the given slave
1639 * element.
1640 */
get_master_dof_indices(EL_DOF_VEC * rvec,const EL_INFO * s_el_info,const FE_SPACE * m_fe_space)1641 const EL_DOF_VEC *get_master_dof_indices(EL_DOF_VEC *rvec,
1642 const EL_INFO *s_el_info,
1643 const FE_SPACE *m_fe_space)
1644 {
1645 FUNCNAME("get_master_dof_indices");
1646 static EL_DOF_VEC *rvec_space;
1647 const BAS_FCTS *m_bfcts = m_fe_space->bas_fcts;
1648 DOF m_dofs[m_bfcts->n_bas_fcts];
1649 const EL *m_el;
1650 const FE_SPACE *m_fe_chain;
1651 int wall;
1652 int orientation;
1653 int type;
1654 int sbf;
1655
1656 DEBUG_TEST_EXIT(s_el_info->fill_flag & FILL_MASTER_INFO,
1657 "slave->master link not set in EL_INFO.\n");
1658
1659 m_el = s_el_info->master.el;
1660 wall = s_el_info->master.opp_vertex;
1661 orientation = s_el_info->master.orientation;
1662 type = s_el_info->master.el_type;
1663
1664 if (INIT_ELEMENT_NEEDED(m_bfcts)) {
1665 EL_INFO m_el_info = { 0 };
1666
1667 fill_master_el_info(&m_el_info, s_el_info, m_bfcts->fill_flags);
1668
1669 if (INIT_ELEMENT(&m_el_info, m_bfcts) == INIT_EL_TAG_NONE) {
1670 return NULL;
1671 }
1672 }
1673
1674 if (rvec == NULL) {
1675 if (rvec_space) {
1676 free_el_dof_vec(rvec_space);
1677 }
1678 rvec = rvec_space = get_el_dof_vec(m_bfcts);
1679 }
1680
1681 GET_DOF_INDICES(m_bfcts, m_el, m_fe_space->admin, m_dofs);
1682 rvec->n_components = m_bfcts->n_trace_bas_fcts[wall];
1683 for (sbf = 0; sbf < rvec->n_components; sbf++) {
1684 rvec->vec[sbf] =
1685 m_dofs[m_bfcts->trace_dof_map[type > 0][orientation < 0][wall][sbf]];
1686 }
1687 CHAIN_FOREACH(m_fe_chain, m_fe_space, const FE_SPACE) {
1688 DOF m_dofs[m_fe_chain->bas_fcts->n_bas_fcts];
1689
1690 rvec = CHAIN_NEXT(rvec, EL_DOF_VEC);
1691
1692 m_bfcts = m_fe_chain->bas_fcts;
1693 GET_DOF_INDICES(m_bfcts, m_el, m_fe_chain->admin, m_dofs);
1694 rvec->n_components = m_bfcts->n_trace_bas_fcts[wall];
1695 for (sbf = 0; sbf < rvec->n_components; sbf++) {
1696 rvec->vec[sbf] =
1697 m_dofs[m_bfcts->trace_dof_map[type > 0][orientation < 0][wall][sbf]];
1698 }
1699 }
1700 rvec = CHAIN_NEXT(rvec, EL_DOF_VEC);
1701
1702 return rvec;
1703 }
1704
1705 /* Compute the boundary classification for all local slave DOFs
1706 * relative to the master mesh.
1707 */
get_master_bound(EL_BNDRY_VEC * rvec,const EL_INFO * s_el_info,const BAS_FCTS * m_bas_fcts)1708 const EL_BNDRY_VEC *get_master_bound(EL_BNDRY_VEC *rvec,
1709 const EL_INFO *s_el_info,
1710 const BAS_FCTS *m_bas_fcts)
1711 {
1712 FUNCNAME("get_master_dof_indices");
1713 static EL_BNDRY_VEC *rvec_space;
1714 BNDRY_FLAGS m_bound[m_bas_fcts->n_bas_fcts_max];
1715 EL_INFO m_el_info = { 0 };
1716 const BAS_FCTS *m_bfcts_chain;
1717 int wall, o, t;
1718 int sbf;
1719
1720 DEBUG_TEST_EXIT(s_el_info->fill_flag & FILL_MASTER_INFO,
1721 "slave->master link not set in EL_INFO.\n");
1722
1723 wall = s_el_info->master.opp_vertex;
1724 o = s_el_info->master.orientation < 0;
1725 t = s_el_info->master.el_type > 0;
1726
1727 fill_master_el_info(&m_el_info, s_el_info, m_bas_fcts->fill_flags|FILL_BOUND);
1728 if (INIT_ELEMENT(&m_el_info, m_bas_fcts) == INIT_EL_TAG_NONE) {
1729 return NULL;
1730 }
1731
1732 if (rvec == NULL) {
1733 if (rvec_space) {
1734 free_el_bndry_vec(rvec_space);
1735 }
1736 rvec = rvec_space = get_el_bndry_vec(m_bas_fcts);
1737 }
1738
1739 GET_BOUND(m_bas_fcts,& m_el_info, m_bound);
1740 rvec->n_components = m_bas_fcts->n_trace_bas_fcts[wall];
1741 for (sbf = 0; sbf < rvec->n_components; sbf++) {
1742 BNDRY_FLAGS_CPY(rvec->vec[sbf],
1743 m_bound[m_bas_fcts->trace_dof_map[t][o][wall][sbf]]);
1744 }
1745 CHAIN_FOREACH(m_bfcts_chain, m_bas_fcts, const BAS_FCTS) {
1746 BNDRY_FLAGS m_bound[m_bfcts_chain->n_bas_fcts];
1747
1748 rvec = CHAIN_NEXT(rvec, EL_BNDRY_VEC);
1749
1750 GET_BOUND(m_bfcts_chain, &m_el_info, m_bound);
1751 rvec->n_components = m_bfcts_chain->n_trace_bas_fcts[wall];
1752 for (sbf = 0; sbf < rvec->n_components; sbf++) {
1753 BNDRY_FLAGS_CPY(rvec->vec[sbf],
1754 m_bound[m_bfcts_chain->trace_dof_map[t][o][wall][sbf]]);
1755 }
1756 }
1757 rvec = CHAIN_NEXT(rvec, EL_BNDRY_VEC);
1758
1759 return rvec;
1760 }
1761
1762 #undef EQ_COPY
1763 #define EQ_COPY(from, to) (to) = (from)
1764
1765 /* A C "template" to compute the trace of a generic DOF_VEC object,
1766 * copyinsn is called via copyinsn(from, to).
1767 */
1768 #define DEFUN_TRACE_DOF_VEC(TYPE, typename, COPYINSN) \
1769 void trace_##typename##_vec(DOF_##TYPE##_VEC *svec, \
1770 const DOF_##TYPE##_VEC *mvec) \
1771 { \
1772 FUNCNAME("trace_"#typename); \
1773 const BAS_FCTS *sbfcts = svec->fe_space->bas_fcts; \
1774 const BAS_FCTS *mbfcts = mvec->fe_space->bas_fcts; \
1775 const DOF_ADMIN *sadmin = svec->fe_space->admin; \
1776 const DOF *s_dofs; \
1777 const EL_DOF_VEC *sm_dofs; \
1778 int sbf; \
1779 \
1780 TEST_EXIT(sbfcts == mbfcts->trace_bas_fcts, \
1781 "svec->fe_space->bas_fcts" \
1782 " != " \
1783 "mvec->fe_space->bas_fcts->trace_bas_fcts!\n"); \
1784 \
1785 TRAVERSE_FIRST(svec->fe_space->mesh, -1, CALL_LEAF_EL|FILL_MASTER_INFO) { \
1786 \
1787 if (INIT_ELEMENT(el_info, sbfcts) == INIT_EL_TAG_NULL) { \
1788 continue; \
1789 } \
1790 \
1791 sm_dofs = get_master_dof_indices(NULL, el_info, mvec->fe_space); \
1792 if (sm_dofs) { \
1793 CHAIN_DO(sm_dofs, EL_DOF_VEC) { \
1794 \
1795 s_dofs = GET_DOF_INDICES(sbfcts, el_info->el, sadmin, NULL)->vec; \
1796 \
1797 for (sbf = 0; sbf < sbfcts->n_bas_fcts; sbf++) { \
1798 COPYINSN(mvec->vec[sm_dofs->vec[sbf]], svec->vec[s_dofs[sbf]]); \
1799 } \
1800 mvec = CHAIN_NEXT(mvec, DOF_##TYPE##_VEC); \
1801 svec = CHAIN_NEXT(svec, DOF_##TYPE##_VEC); \
1802 } CHAIN_WHILE(sm_dofs, EL_DOF_VEC); \
1803 } \
1804 } TRAVERSE_NEXT(); \
1805 }
1806
1807 /* Traces for all kinds of DOF-vectors. :) */
DEFUN_TRACE_DOF_VEC(REAL,dof_real,EQ_COPY)1808 DEFUN_TRACE_DOF_VEC(REAL, dof_real, EQ_COPY)
1809 DEFUN_TRACE_DOF_VEC(REAL_D, dof_real_d, COPY_DOW)
1810 DEFUN_TRACE_DOF_VEC(INT, dof_int, EQ_COPY)
1811 DEFUN_TRACE_DOF_VEC(DOF, dof_dof, EQ_COPY)
1812 DEFUN_TRACE_DOF_VEC(DOF, int_dof, EQ_COPY)
1813 DEFUN_TRACE_DOF_VEC(UCHAR, dof_uchar, EQ_COPY)
1814 DEFUN_TRACE_DOF_VEC(SCHAR, dof_schar, EQ_COPY)
1815 DEFUN_TRACE_DOF_VEC(PTR, dof_ptr, EQ_COPY)
1816
1817 void trace_dof_real_vec_d(DOF_REAL_VEC_D *svec, const DOF_REAL_VEC_D *mvec)
1818 {
1819 FUNCNAME("trace_dof_real_vec_d");
1820 const BAS_FCTS *sbfcts = svec->fe_space->bas_fcts;
1821 const BAS_FCTS *mbfcts = mvec->fe_space->bas_fcts;
1822 const DOF_ADMIN *sadmin = svec->fe_space->admin;
1823 const DOF *s_dofs;
1824 const EL_DOF_VEC *sm_dofs;
1825 int sbf;
1826
1827 TEST_EXIT(sbfcts == mbfcts->trace_bas_fcts,
1828 "svec->fe_space->bas_fcts"
1829 " != "
1830 "mvec->fe_space->bas_fcts->trace_bas_fcts!\n");
1831
1832 TRAVERSE_FIRST(svec->fe_space->mesh, -1, CALL_LEAF_EL|FILL_MASTER_INFO) {
1833
1834 if (INIT_ELEMENT(el_info, sbfcts) == INIT_EL_TAG_NULL) {
1835 continue;
1836 }
1837
1838 sm_dofs = get_master_dof_indices(NULL, el_info, mvec->fe_space);
1839 if (sm_dofs) {
1840 CHAIN_DO(sm_dofs, EL_DOF_VEC) {
1841
1842 s_dofs = GET_DOF_INDICES(sbfcts, el_info->el, sadmin, NULL)->vec;
1843
1844 if (mvec->stride == 1) {
1845 for (sbf = 0; sbf < sbfcts->n_bas_fcts; sbf++) {
1846 svec->vec[s_dofs[sbf]] = mvec->vec[sm_dofs->vec[sbf]];
1847 }
1848 } else {
1849 DOF_REAL_D_VEC *svd = (DOF_REAL_D_VEC *)svec;
1850 DOF_REAL_D_VEC *mvd = (DOF_REAL_D_VEC *)mvec;
1851 for (sbf = 0; sbf < sbfcts->n_bas_fcts; sbf++) {
1852 COPY_DOW(mvd->vec[sm_dofs->vec[sbf]], svd->vec[s_dofs[sbf]]);
1853 }
1854 }
1855 mvec = CHAIN_NEXT(mvec, DOF_REAL_VEC_D);
1856 svec = CHAIN_NEXT(svec, DOF_REAL_VEC_D);
1857 } CHAIN_WHILE(sm_dofs, EL_DOF_VEC);
1858 }
1859 } TRAVERSE_NEXT();
1860 }
1861
1862 /* What follows are some convenience functions to add element matrix
1863 * computed on slave meshes to master matrices. Do not use them. It is
1864 * overkill to derive a slave mesh for the sole purpose of assembling
1865 * boundary integral contributions.
1866 *
1867 * Except when either column of row space live on the trace-mesh
1868 * alone, e.g. for weak slip-conditions.
1869 */
1870
1871 /* Update a DOF_MATRIX living on a master mesh, using an
1872 * EL_MATRIX_INFO object living on a slave mesh.
1873 */
update_master_matrix(DOF_MATRIX * m_dof_matrix,const EL_MATRIX_INFO * s_minfo,MatrixTranspose transpose)1874 void update_master_matrix(DOF_MATRIX *m_dof_matrix,
1875 const EL_MATRIX_INFO *s_minfo,
1876 MatrixTranspose transpose)
1877 {
1878 FUNCNAME("update_master_matrix");
1879 MESH *slave, *row_mesh, *col_mesh;
1880 const FE_SPACE *s_row_fe_space, *s_col_fe_space;
1881 const BAS_FCTS *s_col_bfcts, *m_row_bfcts;
1882 const DOF_ADMIN *s_row_admin;
1883 bool use_get_bound;
1884 EL_DOF_VEC *row_dof, *col_dof;
1885 EL_SCHAR_VEC *bound = NULL;
1886 EL_BNDRY_VEC *bndry_bits = NULL;
1887 FLAGS fill_flag;
1888
1889 TEST_EXIT(s_minfo, "no EL_MATRIX_INFO\n");
1890 TEST_EXIT(s_minfo->el_matrix_fct, "no el_matrix_fct in EL_MATRIX_INFO\n");
1891 TEST_EXIT(m_dof_matrix, "no DOF_MATRIX\n");
1892
1893 slave = s_minfo->row_fe_space->mesh;
1894
1895 BNDRY_FLAGS_CPY(m_dof_matrix->dirichlet_bndry, s_minfo->dirichlet_bndry);
1896
1897 s_col_fe_space = NULL;
1898 if (transpose == NoTranspose) {
1899 s_row_fe_space = s_minfo->row_fe_space;
1900 if (s_minfo->col_fe_space && s_minfo->col_fe_space != s_row_fe_space) {
1901 s_col_fe_space = s_minfo->col_fe_space;
1902 }
1903 } else {
1904 if (s_minfo->col_fe_space &&
1905 s_minfo->col_fe_space != s_minfo->row_fe_space) {
1906 s_row_fe_space = s_minfo->col_fe_space;
1907 s_col_fe_space = s_minfo->row_fe_space;
1908 } else {
1909 s_row_fe_space = s_minfo->col_fe_space;
1910 }
1911 }
1912
1913 s_row_admin = s_row_fe_space->admin;
1914
1915 if (s_col_fe_space) {
1916 s_col_bfcts = s_col_fe_space->bas_fcts;
1917 } else {
1918 s_col_bfcts = NULL;
1919 }
1920
1921 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(m_dof_matrix->dirichlet_bndry);
1922 if (use_get_bound) {
1923 fill_flag = s_minfo->fill_flag|FILL_BOUND;
1924 if (slave->is_periodic && !(s_row_admin->flags & ADM_PERIODIC)) {
1925 fill_flag |= FILL_NON_PERIODIC;
1926 }
1927 } else {
1928 fill_flag = s_minfo->fill_flag;
1929 }
1930
1931 s_minfo->el_matrix_fct(NULL, s_minfo->fill_info);
1932
1933 m_row_bfcts = m_dof_matrix->row_fe_space->bas_fcts;
1934 row_dof = get_el_dof_vec(m_row_bfcts);
1935 if (use_get_bound) {
1936 bound = get_el_schar_vec(m_row_bfcts);
1937 bndry_bits = get_el_bndry_vec(m_row_bfcts);
1938 }
1939
1940 if (s_col_bfcts) {
1941 col_dof = get_el_dof_vec(m_dof_matrix->col_fe_space->bas_fcts);
1942 } else {
1943 col_dof = row_dof;
1944 }
1945
1946 row_mesh = m_dof_matrix->row_fe_space->mesh;
1947 if (m_dof_matrix->col_fe_space != NULL) {
1948 col_mesh = m_dof_matrix->col_fe_space->mesh;
1949 } else {
1950 col_mesh = row_mesh;
1951 }
1952
1953 if (row_mesh != slave && col_mesh != slave) {
1954 TRAVERSE_FIRST(slave, -1, fill_flag) {
1955 const EL_MATRIX *mat;
1956
1957 if ((mat = s_minfo->el_matrix_fct(el_info, s_minfo->fill_info)) == NULL) {
1958 continue;
1959 }
1960
1961 get_master_dof_indices(row_dof, el_info, m_dof_matrix->row_fe_space);
1962 if (s_col_bfcts) {
1963 get_master_dof_indices(col_dof, el_info, m_dof_matrix->col_fe_space);
1964 }
1965 if (use_get_bound) {
1966 get_master_bound(bndry_bits, el_info, m_row_bfcts);
1967 dirichlet_map(bound, bndry_bits, m_dof_matrix->dirichlet_bndry);
1968 }
1969 add_element_matrix(m_dof_matrix, s_minfo->factor, mat, transpose,
1970 row_dof, col_dof, use_get_bound ? bound : NULL);
1971 } TRAVERSE_NEXT();
1972 } else if (row_mesh != slave) {
1973 TRAVERSE_FIRST(slave, -1, fill_flag) {
1974 const EL_MATRIX *mat;
1975
1976 if ((mat = s_minfo->el_matrix_fct(el_info, s_minfo->fill_info)) == NULL) {
1977 continue;
1978 }
1979
1980 get_master_dof_indices(row_dof, el_info, m_dof_matrix->row_fe_space);
1981 if (s_col_bfcts) {
1982 get_dof_indices(col_dof, s_col_fe_space, el_info->el);
1983 }
1984 if (use_get_bound) {
1985 get_master_bound(bndry_bits, el_info, m_row_bfcts);
1986 dirichlet_map(bound, bndry_bits, m_dof_matrix->dirichlet_bndry);
1987 }
1988 add_element_matrix(m_dof_matrix, s_minfo->factor, mat, transpose,
1989 row_dof, col_dof, use_get_bound ? bound : NULL);
1990 } TRAVERSE_NEXT();
1991 } else if (col_mesh != slave) {
1992 TRAVERSE_FIRST(slave, -1, fill_flag) {
1993 const EL_MATRIX *mat;
1994
1995 if ((mat = s_minfo->el_matrix_fct(el_info, s_minfo->fill_info)) == NULL) {
1996 continue;
1997 }
1998
1999 get_dof_indices(row_dof, m_dof_matrix->row_fe_space, el_info->el);
2000 if (s_col_bfcts) {
2001 get_master_dof_indices(col_dof, el_info, m_dof_matrix->col_fe_space);
2002 }
2003 if (use_get_bound) {
2004 get_master_bound(bndry_bits, el_info, m_row_bfcts);
2005 dirichlet_map(bound, bndry_bits, m_dof_matrix->dirichlet_bndry);
2006 }
2007 add_element_matrix(m_dof_matrix, s_minfo->factor, mat, transpose,
2008 row_dof, col_dof, use_get_bound ? bound : NULL);
2009 } TRAVERSE_NEXT();
2010 }
2011
2012 free_el_dof_vec(row_dof);
2013 if (s_col_bfcts) {
2014 free_el_dof_vec(col_dof);
2015 }
2016
2017 if (use_get_bound) {
2018 free_el_schar_vec(bound);
2019 free_el_bndry_vec(bndry_bits);
2020 }
2021 }
2022
update_master_real_vec(DOF_REAL_VEC * m_drv,const EL_VEC_INFO * s_vec_info)2023 void update_master_real_vec(DOF_REAL_VEC *m_drv, const EL_VEC_INFO *s_vec_info)
2024 {
2025 FUNCNAME("update_master_real_vec");
2026 MESH *slave;
2027 const DOF_ADMIN *admin;
2028 bool use_get_bound;
2029 EL_DOF_VEC *dof;
2030 EL_SCHAR_VEC *bound;
2031 const EL_BNDRY_VEC *bndry_bits;
2032 FLAGS fill_flag;
2033
2034 TEST_EXIT(s_vec_info,"no EL_VEC_INFO\n");
2035 TEST_EXIT(s_vec_info->el_vec_fct,"no el_vec_fct in EL_VEC_INFO\n");
2036 TEST_EXIT(m_drv,"no DOF_REAL_VEC\n");
2037
2038 slave = s_vec_info->fe_space->mesh;
2039 admin = s_vec_info->fe_space->admin;
2040
2041 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(s_vec_info->dirichlet_bndry);
2042 if (use_get_bound) {
2043 fill_flag = s_vec_info->fill_flag|FILL_BOUND;
2044 if (slave->is_periodic && !(admin->flags & ADM_PERIODIC)) {
2045 fill_flag |= FILL_NON_PERIODIC;
2046 }
2047 } else {
2048 fill_flag = s_vec_info->fill_flag;
2049 }
2050
2051 s_vec_info->el_vec_fct(NULL, s_vec_info->fill_info);
2052
2053 dof = get_el_dof_vec(m_drv->fe_space->bas_fcts);
2054 bound = get_el_schar_vec(m_drv->fe_space->bas_fcts);
2055
2056 TRAVERSE_FIRST(slave, -1, fill_flag) {
2057 const EL_REAL_VEC *vec;
2058
2059 vec = s_vec_info->el_vec_fct(el_info, s_vec_info->fill_info);
2060 if (vec == NULL) {
2061 continue;
2062 }
2063
2064 get_master_dof_indices(dof, el_info, m_drv->fe_space);
2065 if (use_get_bound) {
2066 bndry_bits = get_master_bound(NULL, el_info, m_drv->fe_space->bas_fcts);
2067 dirichlet_map(bound, bndry_bits, s_vec_info->dirichlet_bndry);
2068 }
2069 add_element_vec(m_drv, s_vec_info->factor, vec, dof,
2070 use_get_bound ? bound : NULL);
2071 } TRAVERSE_NEXT();
2072
2073 free_el_dof_vec(dof);
2074 free_el_schar_vec(bound);
2075 }
2076
update_master_real_d_vec(DOF_REAL_D_VEC * m_drdv,const EL_VEC_D_INFO * s_vec_info)2077 void update_master_real_d_vec(DOF_REAL_D_VEC *m_drdv,
2078 const EL_VEC_D_INFO *s_vec_info)
2079 {
2080 FUNCNAME("update_master_real_d_vec");
2081 MESH *slave;
2082 const BAS_FCTS *bfcts;
2083 const DOF_ADMIN *admin;
2084 bool use_get_bound;
2085 EL_DOF_VEC *dof;
2086 EL_SCHAR_VEC *bound;
2087 EL_BNDRY_VEC *bndry_bits;
2088 FLAGS fill_flag;
2089
2090 TEST_EXIT(s_vec_info,"no EL_VEC_D_INFO\n");
2091 TEST_EXIT(s_vec_info->el_vec_fct,"no el_vec_fct in EL_VEC_D_INFO\n");
2092 TEST_EXIT(m_drdv,"no DOF_REAL_D_VEC\n");
2093
2094 slave = s_vec_info->fe_space->mesh;
2095 bfcts = s_vec_info->fe_space->bas_fcts;
2096 admin = s_vec_info->fe_space->admin;
2097
2098 use_get_bound = !BNDRY_FLAGS_IS_INTERIOR(s_vec_info->dirichlet_bndry);
2099 if (use_get_bound) {
2100 fill_flag = s_vec_info->fill_flag|FILL_BOUND;
2101 if (slave->is_periodic && !(admin->flags & ADM_PERIODIC)) {
2102 fill_flag |= FILL_NON_PERIODIC;
2103 }
2104 } else {
2105 fill_flag = s_vec_info->fill_flag;
2106 }
2107
2108 s_vec_info->el_vec_fct(NULL, s_vec_info->fill_info);
2109
2110 dof = get_el_dof_vec(m_drdv->fe_space->bas_fcts);
2111 bound = get_el_schar_vec(m_drdv->fe_space->bas_fcts);
2112
2113 TRAVERSE_FIRST(slave, -1, fill_flag) {
2114 const EL_REAL_D_VEC *vec;
2115
2116 vec = s_vec_info->el_vec_fct(el_info, s_vec_info->fill_info);
2117 if (vec == NULL) {
2118 continue;
2119 }
2120
2121 get_master_dof_indices(dof, el_info, m_drdv->fe_space);
2122 if (use_get_bound) {
2123 bndry_bits = get_bound(NULL, bfcts, el_info);
2124 dirichlet_map(bound, bndry_bits, s_vec_info->dirichlet_bndry);
2125 }
2126 add_element_d_vec(m_drdv, s_vec_info->factor, vec, dof,
2127 use_get_bound ? bound : NULL);
2128 } TRAVERSE_NEXT();
2129
2130 free_el_dof_vec(dof);
2131 free_el_schar_vec(bound);
2132 }
2133