1 /****************************************************************************/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /* file: HB_precon.c */
7 /* */
8 /* description: hierachical basis preconditioning, including higher */
9 /* order */
10 /* and also the BPX preconditioner */
11 /* */
12 /*--------------------------------------------------------------------------*/
13 /* */
14 /* authors: Alfred Schmidt */
15 /* Zentrum fuer Technomathematik */
16 /* Fachbereich 3 Mathematik/Informatik */
17 /* Universitaet Bremen */
18 /* Bibliothekstr. 2 */
19 /* D-28359 Bremen, Germany */
20 /* */
21 /* Kunibert G. Siebert */
22 /* Institut fuer Mathematik */
23 /* Universitaet Augsburg */
24 /* Universitaetsstr. 14 */
25 /* D-86159 Augsburg, Germany */
26 /* */
27 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
28 /* */
29 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
30 /* */
31 /*--------------------------------------------------------------------------*/
32
33 /*---8<---------------------------------------------------------------------*/
34 /*--- to do: something is wrong in BPX for p > 2 ---*/
35 /*--------------------------------------------------------------------->8---*/
36
37 #ifdef HAVE_CONFIG_H
38 # include "config.h"
39 #endif
40
41 #include "alberta_intern.h"
42 #include "alberta.h"
43
44 typedef struct
45 {
46 PRECON precon;
47 const DOF_MATRIX *matrix;
48 const FE_SPACE *fe_space;
49 const DOF_SCHAR_VEC *bound_dv;
50 BNDRY_FLAGS dirichlet_bndry;
51
52 int info;
53 bool high_degree;
54 int mg_levels;
55 int size;
56 U_CHAR *dof_level;
57 U_CHAR *local_dof;
58 DOF (*dof_parent)[N_VERTICES_MAX];
59 DOF *sort_dof, *sort_dof_invers;
60 int *dofs_per_level;
61 const S_CHAR *bound;
62 REAL (*ipol)[N_VERTICES_MAX]; /* for high_degree interpol. */
63
64 /* BPX section: */
65 REAL *g;
66 REAL diam;
67
68 /* memory stuff */
69 SCRATCH_MEM scratch;
70
71 } HB_DATA;
72
73 static void exit_HB_precon(void *precon_data);
74
75 /****************************************************************************/
76
77 typedef struct hb_traverse_data {
78 int n0_vert, max_level, max_dof_level;
79 const int *n_dof, *n0_dof, *node;
80 int *local_dof_sort;
81 HB_DATA *hb_data;
82 } HB_TRAVERSE_DATA;
83
max_level_fct(const EL_INFO * el_info,void * data)84 static void max_level_fct(const EL_INFO *el_info, void *data)
85 {
86 HB_TRAVERSE_DATA *ud = (HB_TRAVERSE_DATA *)data;
87 int dof, dof0, dof1, level = (int)(el_info->level);
88 int dim = el_info->mesh->dim;
89 EL *el = el_info->el;
90
91 ud->max_level = MAX(level, ud->max_level);
92
93 if (!(IS_LEAF_EL(el))) {
94 dof = el->child[0]->dof[N_VERTICES(dim)-1][ud->n0_vert];
95 dof0 = ud->hb_data->dof_parent[dof][0] = el->dof[0][ud->n0_vert];
96 dof1 = ud->hb_data->dof_parent[dof][1] = el->dof[1][ud->n0_vert];
97 level = 1 + MAX(ud->hb_data->dof_level[dof0], ud->hb_data->dof_level[dof1]);
98 ud->hb_data->dof_level[dof] = level;
99 ud->max_dof_level = MAX(level, ud->max_dof_level);
100 }
101 }
102
103 /* high_degree_fct(): set level of all non-vertex dofs to highest level, */
104 /* and save parent data */
105
high_degree_fct(const EL_INFO * el_info,void * data)106 static void high_degree_fct(const EL_INFO *el_info, void *data)
107 {
108 FUNCNAME("high_degree_fct");
109 HB_TRAVERSE_DATA *ud = (HB_TRAVERSE_DATA *)data;
110 int i, j, k, m, n, nod0, n0, dof;
111 int dim = el_info->mesh->dim;
112 EL *el = el_info->el;
113 DOF vertex_dof[N_VERTICES_MAX];
114
115 DEBUG_TEST_EXIT(IS_LEAF_EL(el), "Non-leaf element???\n");
116
117 for (k=0; k<N_VERTICES(dim); k++)
118 vertex_dof[k] = el->dof[k][ud->n0_vert];
119 m = N_VERTICES(dim);
120
121 if ((n = ud->n_dof[CENTER]) > 0) {
122 nod0 = ud->node[CENTER];
123 n0 = ud->n0_dof[CENTER];
124 for (j=0; j<n; j++) {
125 dof = el->dof[nod0][n0+j];
126 /* MSG("center dof %2d: level %d\n", dof, max_dof_level); */
127 ud->hb_data->dof_level[dof] = ud->max_dof_level;
128 for (k=0; k<N_VERTICES(dim); k++)
129 ud->hb_data->dof_parent[dof][k] = vertex_dof[k];
130 ud->hb_data->local_dof[dof] = ud->local_dof_sort[m++];
131 }
132 }
133
134 if(dim > 1 && (n = ud->n_dof[EDGE]) > 0) {
135 nod0 = ud->node[EDGE];
136 n0 = ud->n0_dof[EDGE];
137 for (i=0; i<N_EDGES(dim); i++)
138 for (j=0; j<n; j++) {
139 dof = el->dof[nod0+i][n0+j];
140 if (!ud->hb_data->local_dof[dof]) {
141 /* MSG("edge dof %2d: level %d\n", dof, max_dof_level); */
142 ud->hb_data->dof_level[dof] = ud->max_dof_level;
143 for (k=0; k<N_VERTICES(dim); k++)
144 ud->hb_data->dof_parent[dof][k] = vertex_dof[k];
145 ud->hb_data->local_dof[dof] = ud->local_dof_sort[m];
146 }
147 m++;
148 }
149
150 }
151
152 #if DIM_OF_WORLD >= 3
153 if (dim == 3 && (n = ud->n_dof[FACE]) > 0) {
154 nod0 = ud->node[FACE];
155 n0 = ud->n0_dof[FACE];
156 for (i=0; i < N_FACES_3D; i++)
157 for (j=0; j<n; j++) {
158 dof = el->dof[nod0+i][n0+j];
159 if (!ud->hb_data->local_dof[dof]) {
160 /* MSG("face dof %2d: level %d\n", dof, max_dof_level); */
161 ud->hb_data->dof_level[dof] = ud->max_dof_level;
162 for (k=0; k < N_VERTICES_3D; k++)
163 ud->hb_data->dof_parent[dof][k] = vertex_dof[k];
164 ud->hb_data->local_dof[dof] = ud->local_dof_sort[m];
165 }
166 m++;
167 }
168 }
169 #endif
170
171 DEBUG_TEST_EXIT(m == ud->hb_data->fe_space->bas_fcts->n_bas_fcts,
172 "m <> n_bas_fcts: %d %d\n",
173 m, ud->hb_data->fe_space->bas_fcts->n_bas_fcts);
174
175 return;
176 }
177
178 /****************************************************************************/
179
180 /* NOTE: This terrible hack is definitely something we could do better... */
181 static DOF *el_hat_dofs[N_VERTICES_MAX+N_EDGES_MAX+N_FACES_MAX+1];
182
183 static MESH dummy_mesh[4] = {{0},{NULL,1},{NULL,2},{NULL,3}};
184
185 static EL el_hat =
186 {
187 {NULL, NULL}, /* children */
188 el_hat_dofs,
189 0 /* ... */
190 };
191
192 #if DIM_MAX == 1
193 static const EL_INFO el_hat_info[2] = {
194 { NULL, /* mesh */
195 {{0.0}, {0.0}}, /* coord */
196 NULL, /* macro el */
197 NULL, NULL, /* el, parent */
198 0, /* fill_flag */
199 0, /* level */
200 { 0 } /*--- and the rest is initialized with 0s ---*/
201 },
202 { dummy_mesh + 1,
203 {{0.0}, {1.0}},
204 NULL,
205 NULL, NULL,
206 FILL_COORDS,
207 0,
208 { 0 } /*--- and the rest is initialized with 0s ---*/
209 }};
210 #elif DIM_MAX == 2
211 static const EL_INFO el_hat_info[3] = {
212 { NULL,
213 {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}},
214 NULL,
215 NULL, NULL,
216 0,
217 0,
218 { 0 } /*--- and the rest is initialized with 0s ---*/
219 },
220 { dummy_mesh + 1,
221 {{0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}},
222 NULL,
223 NULL, NULL,
224 FILL_COORDS,
225 0,
226 { 0 } /*--- and the rest is initialized with 0s ---*/
227 },
228 { dummy_mesh + 2,
229 { {1.0, 0.0}, {0.0, 1.0}, {0.0, 0.0}},
230 NULL,
231 NULL, NULL,
232 FILL_COORDS,
233 0,
234 { 0 } /*--- and the rest is initialized with 0s ---*/
235 }};
236 #elif DIM_MAX >= 3
237 static const EL_INFO el_hat_info[4] = {
238 { NULL,
239 {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}},
240 NULL,
241 NULL, NULL,
242 0,
243 0,
244 { 0 } /*--- and the rest is initialized with 0s ---*/
245 },
246 { dummy_mesh + 1,
247 {{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}},
248 NULL,
249 NULL, NULL,
250 FILL_COORDS,
251 0,
252 { 0 } /*--- and the rest is initialized with 0s ---*/
253 },
254 { dummy_mesh + 2,
255 { {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.0} },
256 NULL,
257 NULL, NULL,
258 FILL_COORDS,
259 0,
260 { 0 } /*--- and the rest is initialized with 0s ---*/
261 },
262 { dummy_mesh + 3,
263 { {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, {0.0, 0.0, 0.0} },
264 NULL,
265 NULL, NULL,
266 FILL_COORDS,
267 0,
268 { 0 } /*--- and the rest is initialized with 0s ---*/
269 }};
270 #endif
271
lambda_0(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)272 static REAL lambda_0(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
273 {
274 return quad->lambda[iq][0];
275 }
lambda_1(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)276 static REAL lambda_1(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
277 {
278 #if DIM_MAX > 0
279 return quad->lambda[iq][1];
280 #else
281 return 0.0;
282 #endif
283 }
lambda_2(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)284 static REAL lambda_2(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
285 {
286 #if DIM_MAX > 1
287 return quad->lambda[iq][2];
288 #else
289 return 0.0;
290 #endif
291 }
lambda_3(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)292 static REAL lambda_3(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
293 {
294 #if DIM_MAX > 2
295 return quad->lambda[iq][3];
296 #else
297 return 0.0;
298 #endif
299 }
300 static REAL (*lambda[N_LAMBDA_LIMIT])(const EL_INFO *el_info,
301 const QUAD *quad, int iq, void *ud) =
302 {
303 lambda_0, lambda_1, lambda_2, lambda_3
304 };
305
306 /****************************************************************************/
307 /* init_HB_precon(matrix, boundary, info); initialize HB_DATA structure */
308 /* fills (void *)(HB_DATA *), if ok, */
309 /* returns NULL, if anything goes wrong. */
310 /****************************************************************************/
311
init_HB_BPX_precon(void * precon_data,int BPX)312 static bool init_HB_BPX_precon(void *precon_data, int BPX)
313 {
314 FUNCNAME("init_HB_BPX_precon");
315 static HB_TRAVERSE_DATA td[1];
316 HB_DATA *data = (HB_DATA *)precon_data;
317 const FE_SPACE *fe_space;
318 const DOF_ADMIN *admin;
319 int i, j, k, m, size, info, dim, hole_count;
320 int *tmp_per_level = NULL;
321 int n_bas_fcts;
322 DOF *tmp_dofs, *tmp_dof_ptr;
323 const DOF *dof_indices;
324 DOF_SCHAR_VEC *dsv;
325
326 if (!data) {
327 ERROR("no precon_data\n");
328 return false;
329 }
330 if (!(data->fe_space)) {
331 ERROR("no precon_data->fe_space\n");
332 return false;
333 }
334
335 info = data->info;
336 fe_space = data->fe_space;
337 if (!fe_space || !(fe_space->admin) || !(fe_space->bas_fcts)) {
338 MSG("no fe_space or admin or bas_fcts.\n");
339 return false;
340 }
341
342 if (fe_space->bas_fcts->n_dof[VERTEX] != 1) {
343 MSG("sorry, only for FE spaces with n_dof[VERTEX]==1.\n");
344 return false;
345 }
346
347 #if 0
348 /* Holes are masked out by the "boundary" mask vector. */
349 TEST_EXIT(fe_space->admin->hole_count == 0,
350 "HB-/BPX-precon does not works with discontiguous index ranges, "
351 "please call dof_compress() first.\n");
352 #endif
353
354 admin = fe_space->admin;
355 size = admin->size_used;
356 dim = fe_space->mesh->dim;
357
358 INIT_OBJECT(fe_space->bas_fcts);
359
360 n_bas_fcts = fe_space->bas_fcts->n_bas_fcts;
361 data->high_degree = (n_bas_fcts > N_VERTICES(dim));
362 if (data->high_degree) {
363 INFO(info, 1, "use high degree version\n");
364 }
365
366 data->mg_levels = 0;
367 data->dofs_per_level = NULL;
368
369 hole_count = 0;
370 dsv = get_dof_schar_vec_skel("HB/BPX bound", data->fe_space, data->scratch);
371 data->bound = dsv->vec = SCRATCH_MEM_ALLOC(data->scratch, dsv->size, S_CHAR);
372 if (data->bound_dv) {
373 for (i = 0; i < dsv->size; i++) {
374 if (data->matrix->matrix_row[i] == NULL) {
375 dsv->vec[i] = DIRICHLET;
376 ++hole_count;
377 } else {
378 dsv->vec[i] = data->bound_dv->vec[i];
379 }
380 }
381 } else if (!BNDRY_FLAGS_IS_INTERIOR(data->matrix->dirichlet_bndry)) {
382 dirichlet_bound(NULL, NULL, dsv, data->dirichlet_bndry, NULL);
383 for (i = 0; i < dsv->size; i++) {
384 if (data->matrix->matrix_row[i] == NULL) {
385 dsv->vec[i] = DIRICHLET;
386 ++hole_count;
387 }
388 }
389 } else {
390 for (i = 0; i < dsv->size; i++) {
391 if (data->matrix->matrix_row[i] == NULL) {
392 dsv->vec[i] = DIRICHLET;
393 ++hole_count;
394 } else {
395 dsv->vec[i] = INTERIOR;
396 }
397 }
398 }
399 if (hole_count == 0 &&
400 data->bound_dv == NULL &&
401 BNDRY_FLAGS_IS_INTERIOR(data->matrix->dirichlet_bndry)) {
402 data->bound = NULL; /* No DOF-mask needed. */
403 }
404
405 data->ipol = NULL;
406
407 data->dof_level = SCRATCH_MEM_ALLOC(data->scratch, 2*size, U_CHAR);
408 data->local_dof = data->dof_level + size;
409
410 data->sort_dof = SCRATCH_MEM_ALLOC(data->scratch, (N_VERTICES_MAX+2)*size, DOF);
411 data->sort_dof_invers = data->sort_dof + size;
412 data->dof_parent = (DOF (*)[N_VERTICES_MAX])(data->sort_dof_invers + size);
413 data->size = size;
414
415 FOR_ALL_DOFS(admin,
416 data->dof_level[dof] = 0;
417 for (j=0; j < N_VERTICES_MAX;j++) data->dof_parent[dof][j] = -1;
418 data->local_dof[dof] = 0;
419 );
420
421 td->n0_vert = admin->n0_dof[VERTEX];
422 td->max_level = td->max_dof_level = 0;
423 td->hb_data = data;
424 mesh_traverse(fe_space->mesh, -1, CALL_EVERY_EL_PREORDER, max_level_fct, td);
425 data->mg_levels = (td->max_level + dim - 1) / dim + 1;
426
427 TEST_EXIT(data->mg_levels == (td->max_dof_level+1),
428 "mg_levels %d != max_dof_level %d + 1\n",
429 data->mg_levels, td->max_dof_level);
430
431 if (data->high_degree) { /* add level for fine-grid non-vertex DOFs */
432 DEF_EL_VEC_VAR(REAL, one_ipol, n_bas_fcts, n_bas_fcts, false);
433
434 data->mg_levels++;
435 data->ipol = (REAL (*)[N_VERTICES_MAX])
436 (SCRATCH_MEM_ALLOC(data->scratch, N_VERTICES_MAX*n_bas_fcts, REAL));
437
438 td->local_dof_sort = SCRATCH_MEM_ALLOC(data->scratch, n_bas_fcts, int);
439 for (j=0; j<n_bas_fcts; j++) {
440 td->local_dof_sort[j] = j; /* ???? */
441 }
442
443 tmp_dof_ptr = tmp_dofs = SCRATCH_MEM_ALLOC(data->scratch, fe_space->mesh->n_dof_el, DOF);
444 m = 0;
445 for (i=0; i<N_VERTICES(dim); i++) {
446 el_hat_dofs[fe_space->mesh->node[VERTEX]+i] = tmp_dof_ptr;
447 for (j=0; j<admin->n_dof[VERTEX]; j++) {
448 tmp_dof_ptr[admin->n0_dof[VERTEX]+j] = m++;
449 }
450 /* If n0_dof[VERTEX] != 0, then we should also fill
451 * tmp_dof_ptr[0+j], otherwise sort_face_indices_3d() will
452 * become confused.
453 */
454 if (admin->n0_dof[VERTEX] != 0) {
455 for (j=0; j < admin->n_dof[VERTEX]; j++) {
456 tmp_dof_ptr[0+j] = tmp_dof_ptr[admin->n0_dof[VERTEX]+j];
457 }
458 }
459 tmp_dof_ptr += fe_space->mesh->n_dof[VERTEX];
460 }
461
462 if (fe_space->mesh->n_dof[CENTER]) {
463 el_hat_dofs[fe_space->mesh->node[CENTER]] = tmp_dof_ptr;
464 for (j=0; j<admin->n_dof[CENTER]; j++) {
465 tmp_dof_ptr[admin->n0_dof[CENTER]+j] = m++;
466 }
467 }
468
469 if(dim > 1 && fe_space->mesh->n_dof[EDGE])
470 for (i=0; i < N_EDGES(dim); i++) {
471 el_hat_dofs[fe_space->mesh->node[EDGE]+i] = tmp_dof_ptr;
472 for (j=0; j<admin->n_dof[EDGE]; j++) {
473 tmp_dof_ptr[admin->n0_dof[EDGE]+j] = m++;
474 }
475 tmp_dof_ptr += fe_space->mesh->n_dof[EDGE];
476 }
477
478 if (dim == 3 && fe_space->mesh->n_dof[FACE])
479 for (i=0; i<N_FACES_3D; i++) {
480 el_hat_dofs[fe_space->mesh->node[FACE]+i] = tmp_dof_ptr;
481 for (j=0; j<admin->n_dof[FACE]; j++) {
482 tmp_dof_ptr[admin->n0_dof[FACE]+j] = m++;
483 }
484 tmp_dof_ptr += fe_space->mesh->n_dof[FACE];
485 }
486
487 TEST_EXIT(m==n_bas_fcts,"m != n_bas_fcts: %d %d\n", m, n_bas_fcts);
488
489 dof_indices = GET_DOF_INDICES(fe_space->bas_fcts, &el_hat, admin, NULL)->vec;
490 for (i = 0; i<fe_space->bas_fcts->n_bas_fcts; i++) {
491 td->local_dof_sort[dof_indices[i]] = i;
492 }
493
494 #if 0
495 print_int_vec("dof_indices ", dof_indices, n_bas_fcts);
496 print_int_vec("local_dof_sort", td->local_dof_sort, n_bas_fcts);
497 #endif
498
499 for (i=0; i < N_VERTICES(dim); i++) {
500 INTERPOL(fe_space->bas_fcts, one_ipol,
501 &(el_hat_info[dim]), -1, 0, NULL, lambda[i], NULL);
502 for (j=0; j<n_bas_fcts; j++)
503 data->ipol[j][i] = one_ipol->vec[j];
504 }
505
506 td->max_dof_level++;
507 td->n0_dof = admin->n0_dof;
508 td->n_dof = admin->n_dof;
509 td->node = fe_space->mesh->node;
510
511 mesh_traverse(fe_space->mesh, -1, CALL_LEAF_EL, high_degree_fct, td);
512
513 if (info > 3)
514 {
515 for (i=0; i<n_bas_fcts; i++)
516 {
517 MSG("ipol[%2d]", i); PRINT_REAL_VEC("", data->ipol[i],N_VERTICES(dim));
518 }
519 }
520
521 #if 0
522 FOR_ALL_DOFS(admin,
523 if (data->local_dof[dof])
524 MSG("dof %3d: local_dof=%2d\n", dof, data->local_dof[dof]);
525 );
526 #endif
527 }
528
529 if (data->mg_levels < 2) {
530 /* exit_HB_precon(data); */
531 return false;
532 }
533
534 data->dofs_per_level = SCRATCH_MEM_ALLOC(data->scratch, data->mg_levels, int);
535 tmp_per_level = SCRATCH_MEM_ALLOC(data->scratch, data->mg_levels, int);
536
537 for (i = 0; i < data->mg_levels; i++) data->dofs_per_level[i] = 0;
538 FOR_ALL_DOFS(admin, data->dofs_per_level[data->dof_level[dof]]++);
539 if (info > 3) {
540 MSG("dofs_per_level:");
541 for (i = 0; i < data->mg_levels; i++)
542 print_msg(" %d", data->dofs_per_level[i]);
543 print_msg("\n");
544 }
545 for (i = 1; i < data->mg_levels; i++) {
546 tmp_per_level[i] = data->dofs_per_level[i-1];
547 data->dofs_per_level[i] += data->dofs_per_level[i-1];
548 }
549 tmp_per_level[0] = 0; /* pointers for filling the sort vectors */
550
551 if (info > 3) {
552 MSG("dofs_per_level accumulated:");
553 for (i = 0; i < data->mg_levels; i++)
554 print_msg(" %d", data->dofs_per_level[i]);
555 print_msg("\n");
556 }
557
558 #if 0
559 for (i = 0; i < data->dofs_per_level[data->mg_levels-1]; i++)
560 MSG("dof_parent[%3d] = (%3d,%3d,%3d), lev=%2d (%2d,%2d,%3d)\n",
561 i,
562 data->dof_parent[i][0], data->dof_parent[i][1], data->dof_parent[i][2],
563 data->dof_level[i],
564 data->dof_level[data->dof_parent[i][0]],
565 data->dof_level[data->dof_parent[i][1]],
566 data->dof_level[data->dof_parent[i][2]]);
567 #endif
568
569 /* build sort_dof[] and sort_dof_invers[] vectors */
570
571 FOR_ALL_DOFS(admin,
572 j = data->dof_level[dof];
573 k = tmp_per_level[j]++;
574 data->sort_dof[k] = dof;
575 data->sort_dof_invers[dof] = k;
576 );
577
578 #if 0
579 for (i = 0; i < data->dofs_per_level[data->mg_levels-1]; i++)
580 {
581 j = data->sort_dof[i];
582 MSG("sort[%3d]: dof=%3d, lev=%2d; invers[%3d]=%3d\n",
583 i, j, data->dof_level[j], j, data->sort_dof_invers[j]);
584 }
585 /* WAIT; */
586 #endif
587
588 if (BPX) {
589 data->g = SCRATCH_MEM_CALLOC(data->scratch, DIM_OF_WORLD*data->size, REAL);
590 data->diam = fe_space->mesh->diam[0];
591 for (i=1; i<DIM_OF_WORLD; i++)
592 data->diam = MAX(data->diam, fe_space->mesh->diam[i]);
593 } else {
594 data->g = NULL;
595 data->diam = 0;
596 }
597
598 return true;
599 }
600
601 /****************************************************************************/
602 /****************************************************************************/
603
exit_HB_BPX_precon(void * vdata,int BPX)604 static void exit_HB_BPX_precon(void *vdata, int BPX)
605 {
606 FUNCNAME("exit_HB_BPX_precon");
607 HB_DATA *data = (HB_DATA *)vdata;
608
609 if (!data) {
610 MSG("no data ???\n");
611 return;
612 }
613
614 SCRATCH_MEM_ZAP(data->scratch);
615 }
616
617 /****************************************************************************/
618
init_HB_precon(void * precon_data)619 static bool init_HB_precon(void *precon_data)
620 {
621 return init_HB_BPX_precon(precon_data, 0);
622 }
623
exit_HB_precon(void * precon_data)624 static void exit_HB_precon(void *precon_data)
625 {
626 exit_HB_BPX_precon(precon_data, 0);
627 }
628
629 /****************************************************************************/
630
HB_precon_s(void * vdata,int n,REAL * r)631 static void HB_precon_s(void *vdata, int n, REAL *r)
632 {
633 FUNCNAME("HB_precon_s");
634 int i, idof, level, last, jdof, k, level1, dim;
635 HB_DATA *data;
636
637 data = (HB_DATA *)vdata;
638 if (!data) {
639 MSG("no data ???\n");
640 return;
641 }
642 dim = data->fe_space->mesh->dim;
643
644 if (n > data->size) {
645 MSG("n > data->size ???\n");
646 return;
647 }
648
649 if (data->mg_levels < 2) /* nothing to do */
650 {
651 return;
652 }
653
654 /* transposed basis transformation (run over refined levels only) */
655
656 if (data->high_degree) {
657 last = data->dofs_per_level[data->mg_levels - 1];
658 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
659 idof = data->sort_dof[i];
660 jdof = data->local_dof[i];
661 if (data->bound) {
662 for (k=0; k<N_VERTICES(dim); k++) {
663 if (data->bound[data->dof_parent[idof][k]] <= INTERIOR)
664 r[data->dof_parent[idof][k]] += data->ipol[jdof][k] * r[idof];
665 }
666 }
667 else {
668 for (k=0; k<N_VERTICES(dim); k++) {
669 r[data->dof_parent[idof][k]] += data->ipol[jdof][k] * r[idof];
670 }
671 }
672 }
673 level1 = data->mg_levels - 2;
674 }
675 else {
676 level1 = data->mg_levels - 1;
677 }
678
679 for (level = level1; level > 0; level--)
680 {
681 last = data->dofs_per_level[level];
682
683 for (i = data->dofs_per_level[level-1]; i < last; i++)
684 {
685 idof = data->sort_dof[i];
686 if (data->bound) {
687 if (data->bound[data->dof_parent[idof][0]] <= INTERIOR)
688 r[data->dof_parent[idof][0]] += 0.5 * r[idof];
689 if (data->bound[data->dof_parent[idof][1]] <= INTERIOR)
690 r[data->dof_parent[idof][1]] += 0.5 * r[idof];
691 }
692 else {
693 r[data->dof_parent[idof][0]] += 0.5 * r[idof];
694 r[data->dof_parent[idof][1]] += 0.5 * r[idof];
695 }
696 }
697 }
698
699 /* basis transformation (run over refined levels only) */
700
701 for (level = 1; level <= level1; level++)
702 {
703 last = data->dofs_per_level[level];
704
705 for (i = data->dofs_per_level[level-1]; i < last; i++)
706 {
707 idof = data->sort_dof[i];
708 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
709 r[idof] += 0.5 * (r[data->dof_parent[idof][0]]
710 + r[data->dof_parent[idof][1]]);
711 }
712 }
713 }
714
715 if (data->high_degree) {
716 last = data->dofs_per_level[data->mg_levels - 1];
717 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
718 idof = data->sort_dof[i];
719 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
720 jdof = data->local_dof[i];
721 for (k = 0; k < N_VERTICES(dim); k++)
722 r[idof] += data->ipol[jdof][k] * r[data->dof_parent[idof][k]];
723 }
724 }
725 }
726
727 return;
728 }
729
730 /****************************************************************************/
731
HB_precon_d(void * vdata,int n,REAL * r)732 static void HB_precon_d(void *vdata, int n, REAL *r)
733 {
734 FUNCNAME("HB_precon_d");
735 int i, idof, level, last, jdof, k, level1, dim;
736 HB_DATA *data;
737 REAL_D *rd = (REAL_D *)r;
738
739 data = (HB_DATA *)vdata;
740 if (!data) {
741 MSG("no data ???\n");
742 return;
743 }
744 dim = data->fe_space->mesh->dim;
745
746 if (n > DIM_OF_WORLD*data->size) {
747 MSG("n > DIM_OF_WORLD*data->size ???\n");
748 return;
749 }
750
751 if (data->mg_levels < 2) /* nothing to do */
752 {
753 return;
754 }
755
756 /* transposed basis transformation (run over refined levels only) */
757
758 if (data->high_degree) {
759 last = data->dofs_per_level[data->mg_levels - 1];
760 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
761 idof = data->sort_dof[i];
762 jdof = data->local_dof[i];
763 if (data->bound) {
764 for (k=0; k<N_VERTICES(dim); k++) {
765 if (data->bound[data->dof_parent[idof][k]] <= INTERIOR)
766 AXPY_DOW(data->ipol[jdof][k], rd[idof], rd[data->dof_parent[idof][k]]);
767 }
768 }
769 else {
770 for (k=0; k<N_VERTICES(dim); k++) {
771 AXPY_DOW(data->ipol[jdof][k], rd[idof],
772 rd[data->dof_parent[idof][k]]);
773 }
774 }
775 }
776 level1 = data->mg_levels - 2;
777 }
778 else {
779 level1 = data->mg_levels - 1;
780 }
781
782
783 for (level = level1; level > 0; level--)
784 {
785 last = data->dofs_per_level[level];
786
787 for (i = data->dofs_per_level[level-1]; i < last; i++)
788 {
789 idof = data->sort_dof[i];
790 if (data->bound) {
791 if (data->bound[data->dof_parent[idof][0]] <= INTERIOR)
792 AXPY_DOW(0.5, rd[idof], rd[data->dof_parent[idof][0]]);
793 if (data->bound[data->dof_parent[idof][1]] <= INTERIOR)
794 AXPY_DOW(0.5, rd[idof], rd[data->dof_parent[idof][1]]);
795 }
796 else {
797 AXPY_DOW(0.5, rd[idof], rd[data->dof_parent[idof][0]]);
798 AXPY_DOW(0.5, rd[idof], rd[data->dof_parent[idof][1]]);
799 }
800 }
801 }
802
803 /* basis transformation (run over refined levels only) */
804
805 for (level = 1; level <= level1; level++)
806 {
807 last = data->dofs_per_level[level];
808
809 for (i = data->dofs_per_level[level-1]; i < last; i++)
810 {
811 idof = data->sort_dof[i];
812 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
813 AXPBYP_DOW(0.5, rd[data->dof_parent[idof][0]],
814 0.5, rd[data->dof_parent[idof][1]],
815 rd[idof]);
816 }
817 }
818 }
819
820 if (data->high_degree) {
821 last = data->dofs_per_level[data->mg_levels - 1];
822 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
823 idof = data->sort_dof[i];
824 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
825 jdof = data->local_dof[i];
826 for (k = 0; k < N_VERTICES(dim); k++)
827 AXPY_DOW(data->ipol[jdof][k], rd[data->dof_parent[idof][k]],
828 rd[idof]);
829 }
830 }
831 }
832
833 return;
834 }
835
836 /****************************************************************************/
837 /****************************************************************************/
838
get_HB_precon_s(const DOF_MATRIX * matrix,const DOF_SCHAR_VEC * bound,int info)839 static const PRECON *get_HB_precon_s(const DOF_MATRIX *matrix,
840 const DOF_SCHAR_VEC *bound,
841 int info)
842 {
843 PRECON *precon;
844 HB_DATA *data;
845 SCRATCH_MEM scratch;
846 const FE_SPACE *fe_space = matrix->row_fe_space;
847
848 if (bound && !FE_SPACE_EQ_P(bound->fe_space, fe_space)) {
849 ERROR("different fe spaces ?\n");
850 return NULL;
851 }
852
853 #if 0
854 /* this stuff does not work with DOF-space holes */
855 dof_compress(fe_space->mesh);
856 #endif
857
858 SCRATCH_MEM_INIT(scratch);
859 data = SCRATCH_MEM_CALLOC(scratch, 1, HB_DATA);
860 SCRATCH_MEM_CPY(data->scratch, scratch);
861
862 precon = &data->precon;
863
864 data->matrix = matrix;
865 data->fe_space = fe_space;
866 data->bound_dv = bound;
867 BNDRY_FLAGS_CPY(data->dirichlet_bndry, data->matrix->dirichlet_bndry);
868
869 data->info = info;
870
871 precon->precon_data = data;
872 precon->init_precon = init_HB_precon;
873 precon->precon = HB_precon_s;
874 precon->exit_precon = exit_HB_precon;
875
876 return precon;
877 }
878
get_HB_precon_d(const DOF_MATRIX * matrix,const DOF_SCHAR_VEC * bound,int info)879 static const PRECON *get_HB_precon_d(const DOF_MATRIX *matrix,
880 const DOF_SCHAR_VEC *bound,
881 int info)
882 {
883 PRECON *precon;
884 HB_DATA *data;
885 SCRATCH_MEM scratch;
886 const FE_SPACE *fe_space = matrix->row_fe_space;
887
888 if (bound && !FE_SPACE_EQ_P(bound->fe_space, fe_space)) {
889 ERROR("different fe spaces ?\n");
890 return NULL;
891 }
892
893 #if 0
894 /* this stuff does not work with DOF-space holes */
895 dof_compress(fe_space->mesh);
896 #endif
897
898 SCRATCH_MEM_INIT(scratch);
899
900 data = SCRATCH_MEM_CALLOC(scratch, 1, HB_DATA);
901 SCRATCH_MEM_CPY(data->scratch, scratch);
902
903 precon = &data->precon;
904
905 data->matrix = matrix;
906 data->fe_space = fe_space;
907 data->bound_dv = bound;
908 BNDRY_FLAGS_CPY(data->dirichlet_bndry, matrix->dirichlet_bndry);
909
910 data->info = info;
911
912 precon->precon_data = data;
913 precon->init_precon = init_HB_precon;
914 precon->precon = HB_precon_d;
915 precon->exit_precon = exit_HB_precon;
916
917 return precon;
918 }
919
920 /****************************************************************************/
921 /****************************************************************************/
922 /** BPX preconditioner **/
923 /****************************************************************************/
924 /****************************************************************************/
925
init_BPX_precon(void * precon_data)926 static bool init_BPX_precon(void *precon_data)
927 {
928 return init_HB_BPX_precon(precon_data, 1);
929 }
930
exit_BPX_precon(void * precon_data)931 static void exit_BPX_precon(void *precon_data)
932 {
933 exit_HB_BPX_precon(precon_data, 1);
934 }
935 /****************************************************************************/
936
BPX_precon_s(void * vdata,int n,REAL * h)937 static void BPX_precon_s(void *vdata, int n, REAL *h)
938 {
939 FUNCNAME("BPX_precon_s");
940 int i, idof, level, first = 0, last, jdof, k, level1, *dof_parent, dim;
941 HB_DATA *data;
942 REAL *g, *ipol;
943
944 data = (HB_DATA *)vdata;
945 if (!data) {
946 MSG("no data ???\n");
947 return;
948 }
949 dim = data->fe_space->mesh->dim;
950
951 if (n > data->size) {
952 MSG("n > data->size ???\n");
953 return;
954 }
955
956 if (data->mg_levels < 2) /* nothing to do */
957 {
958 return;
959 }
960
961 g = data->g;
962 DEBUG_TEST_EXIT(g, "no g vec in HB_DATA\n");
963
964 /* copy h to g */
965 for (i=0; i<data->size; i++) g[i] = h[i];
966
967 /* Loop over all non-macro levels */
968
969 if (data->high_degree) {
970 last = data->dofs_per_level[data->mg_levels - 1];
971
972 /* inverse basis transformation on high degree level for h */
973
974 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
975 idof = data->sort_dof[i];
976 dof_parent = data->dof_parent[idof];
977 jdof = data->local_dof[i];
978 ipol = data->ipol[jdof];
979 if (data->bound) {
980 for (k=0; k<N_VERTICES(dim); k++) {
981 if (data->bound[data->dof_parent[idof][k]] <= INTERIOR)
982 h[idof] -= ipol[k] * h[dof_parent[k]];
983 }
984 }
985 else {
986 for (k=0; k<N_VERTICES(dim); k++) {
987 h[idof] -= ipol[k] * h[dof_parent[k]];
988 }
989 }
990 }
991
992 /* transposed basis transformation on high degree level for g */
993
994 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
995 idof = data->sort_dof[i];
996 dof_parent = data->dof_parent[idof];
997 jdof = data->local_dof[i];
998 ipol = data->ipol[jdof];
999 if (data->bound) {
1000 for (k=0; k<N_VERTICES(dim); k++) {
1001 if (data->bound[dof_parent[k]] <= INTERIOR)
1002 g[dof_parent[k]] += ipol[k] * g[idof];
1003 }
1004 }
1005 else {
1006 for (k=0; k<N_VERTICES(dim); k++) {
1007 g[dof_parent[k]] += ipol[k] * g[idof];
1008 }
1009 }
1010 }
1011
1012 /* add up to degree1 level */
1013
1014 last = data->dofs_per_level[data->mg_levels - 2];
1015 if (data->bound) {
1016 for (i = 0; i < last; i++) {
1017 idof = data->sort_dof[i];
1018 if (data->bound[idof] <= 0)
1019 h[idof] += g[idof];
1020 }
1021 } else
1022 {
1023 for (i = first; i < last; i++)
1024 {
1025 idof = data->sort_dof[i];
1026 h[idof] += g[idof];
1027 }
1028 }
1029
1030 level1 = data->mg_levels - 2;
1031 } else {
1032 level1 = data->mg_levels - 1;
1033 }
1034
1035 for (level = level1; level > 0; level--)
1036 {
1037 last = data->dofs_per_level[level];
1038
1039 /* inverse basis transformation on level k (only) for h */
1040 for (i = data->dofs_per_level[level-1]; i < last; i++) {
1041 idof = data->sort_dof[i];
1042 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
1043 h[idof] -= 0.5 * (h[data->dof_parent[idof][0]]
1044 + h[data->dof_parent[idof][1]]);
1045 }
1046 }
1047
1048 /* transposed basis transformation on level k-1 for g */
1049 for (i = data->dofs_per_level[level-1]; i < last; i++)
1050 {
1051 idof = data->sort_dof[i];
1052 if (data->bound) {
1053 if (data->bound[data->dof_parent[idof][0]] <= INTERIOR)
1054 g[data->dof_parent[idof][0]] += 0.5 * g[idof];
1055 if (data->bound[data->dof_parent[idof][1]] <= INTERIOR)
1056 g[data->dof_parent[idof][1]] += 0.5 * g[idof];
1057 }
1058 else {
1059 g[data->dof_parent[idof][0]] += 0.5 * g[idof];
1060 g[data->dof_parent[idof][1]] += 0.5 * g[idof];
1061 }
1062 }
1063
1064 /* add up to level k-1 */
1065
1066
1067 last = data->dofs_per_level[level-1];
1068 if (data->bound) {
1069 for (i = 0; i < last; i++) {
1070 idof = data->sort_dof[i];
1071 if (data->bound[idof] <= 0)
1072 h[idof] += g[idof];
1073 }
1074 } else
1075 {
1076 for (i = first; i < last; i++)
1077 {
1078 idof = data->sort_dof[i];
1079 h[idof] += g[idof];
1080 }
1081 }
1082 }
1083
1084 /* basis transformation (all refined levels) */
1085
1086 for (level = 1; level <= level1; level++)
1087 {
1088 last = data->dofs_per_level[level];
1089
1090 for (i = data->dofs_per_level[level-1]; i < last; i++)
1091 {
1092 idof = data->sort_dof[i];
1093 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
1094 h[idof] += 0.5 * (h[data->dof_parent[idof][0]]
1095 + h[data->dof_parent[idof][1]]);
1096 }
1097 }
1098 }
1099
1100 if (data->high_degree) {
1101 last = data->dofs_per_level[data->mg_levels - 1];
1102 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
1103 idof = data->sort_dof[i];
1104 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
1105 dof_parent = data->dof_parent[idof];
1106 jdof = data->local_dof[i];
1107 ipol = data->ipol[jdof];
1108 for (k = 0; k < N_VERTICES(dim); k++)
1109 h[idof] += ipol[k] * h[dof_parent[k]];
1110 }
1111 }
1112 }
1113
1114 return;
1115 }
1116
1117 /****************************************************************************/
BPX_precon_d(void * vdata,int n,REAL * hh)1118 static void BPX_precon_d(void *vdata, int n, REAL *hh)
1119 {
1120 FUNCNAME("BPX_precon_d");
1121 int i, k, idof, level, first = 0, last, jdof, level1, dim;
1122 HB_DATA *data;
1123 REAL_D *g;
1124 REAL_D *h = (REAL_D *)hh;
1125
1126 data = (HB_DATA *)vdata;
1127 if (!data) {
1128 MSG("no data ???\n");
1129 return;
1130 }
1131 dim = data->fe_space->mesh->dim;
1132
1133 if (n > data->size*DIM_OF_WORLD) {
1134 MSG("n > data->size*DIM_OF_WORLD ???\n");
1135 return;
1136 }
1137
1138 if (data->mg_levels < 2) /* nothing to do */
1139 {
1140 return;
1141 }
1142
1143 g = (REAL_D *)data->g;
1144 DEBUG_TEST_EXIT(g, "no g vec in HB_DATA\n");
1145
1146 /* copy h to g */
1147 for (i=0; i<data->size; i++)
1148 COPY_DOW(h[i], g[i]);
1149
1150 /* Loop over all non-macro levels */
1151
1152 if (data->high_degree) {
1153 last = data->dofs_per_level[data->mg_levels - 1];
1154
1155 /* inverse basis transformation on high degree level for h */
1156
1157 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
1158 idof = data->sort_dof[i];
1159 jdof = data->local_dof[i];
1160 if (data->bound) {
1161 for (k=0; k<N_VERTICES(dim); k++) {
1162 if (data->bound[data->dof_parent[idof][k]] <= INTERIOR)
1163 AXPY_DOW(-data->ipol[jdof][k], h[data->dof_parent[idof][k]],
1164 h[idof]);
1165 }
1166 }
1167 else {
1168 for (k=0; k<N_VERTICES(dim); k++) {
1169 AXPY_DOW(-data->ipol[jdof][k], h[data->dof_parent[idof][k]], h[idof]);
1170 }
1171 }
1172 }
1173
1174 /* transposed basis transformation on high degree level for g */
1175
1176 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
1177 idof = data->sort_dof[i];
1178 jdof = data->local_dof[i];
1179 if (data->bound) {
1180 for (k=0; k<N_VERTICES(dim); k++) {
1181 if (data->bound[data->dof_parent[idof][k]] <= INTERIOR)
1182 AXPY_DOW(data->ipol[jdof][k], g[idof],
1183 g[data->dof_parent[idof][k]]);
1184 }
1185 }
1186 else {
1187 for (k=0; k<N_VERTICES(dim); k++) {
1188 AXPY_DOW(data->ipol[jdof][k], g[idof], g[data->dof_parent[idof][k]]);
1189 }
1190 }
1191 }
1192
1193 /* add up to degree1 level */
1194
1195 last = data->dofs_per_level[data->mg_levels - 2];
1196 if (data->bound) {
1197 for (i = 0; i < last; i++) {
1198 idof = data->sort_dof[i];
1199 if (data->bound[idof] <= 0)
1200 AXPY_DOW(1.0, g[idof], h[idof]);
1201 }
1202 } else {
1203 for (i = first; i < last; i++) {
1204 idof = data->sort_dof[i];
1205 AXPY_DOW(1.0, g[idof], h[idof]);
1206 }
1207 }
1208
1209 level1 = data->mg_levels - 2;
1210 } else {
1211 level1 = data->mg_levels - 1;
1212 }
1213
1214 for (level = level1; level > 0; level--)
1215 {
1216 last = data->dofs_per_level[level];
1217
1218 /* inverse basis transformation on level k (only) for h */
1219 for (i = data->dofs_per_level[level-1]; i < last; i++) {
1220 idof = data->sort_dof[i];
1221 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
1222 AXPBYP_DOW(0.5, h[data->dof_parent[idof][0]],
1223 0.5, h[data->dof_parent[idof][1]],
1224 h[idof]);
1225 }
1226 }
1227
1228 /* transposed basis transformation on level k-1 for g */
1229 for (i = data->dofs_per_level[level-1]; i < last; i++)
1230 {
1231 idof = data->sort_dof[i];
1232 if (data->bound) {
1233 if (data->bound[data->dof_parent[idof][0]] <= INTERIOR)
1234 AXPY_DOW(0.5, g[idof], g[data->dof_parent[idof][0]]);
1235 if (data->bound[data->dof_parent[idof][1]] <= INTERIOR)
1236 AXPY_DOW(0.5, g[idof], g[data->dof_parent[idof][1]]);
1237 }
1238 else {
1239 AXPY_DOW(0.5, g[idof], g[data->dof_parent[idof][0]]);
1240 AXPY_DOW(0.5, g[idof], g[data->dof_parent[idof][1]]);
1241 }
1242 }
1243
1244 /* add up to level k-1 */
1245
1246 last = data->dofs_per_level[level-1];
1247 if (data->bound) {
1248 for (i = 0; i < last; i++) {
1249 idof = data->sort_dof[i];
1250 if (data->bound[idof] <= 0)
1251 AXPY_DOW(1.0, g[idof], h[idof]);
1252 }
1253 } else {
1254 for (i = first; i < last; i++) {
1255 idof = data->sort_dof[i];
1256 AXPY_DOW(1.0, g[idof], h[idof]);
1257 }
1258 }
1259
1260 }
1261
1262 /* basis transformation (all refined levels) */
1263
1264 for (level = 1; level <= level1; level++)
1265 {
1266 last = data->dofs_per_level[level];
1267
1268 for (i = data->dofs_per_level[level-1]; i < last; i++)
1269 {
1270 idof = data->sort_dof[i];
1271 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
1272 AXPBYP_DOW(0.5, h[data->dof_parent[idof][0]],
1273 0.5, h[data->dof_parent[idof][1]],
1274 h[idof]);
1275 }
1276 }
1277 }
1278
1279 if (data->high_degree) {
1280 last = data->dofs_per_level[data->mg_levels - 1];
1281 for (i = data->dofs_per_level[data->mg_levels - 2]; i < last; i++) {
1282 idof = data->sort_dof[i];
1283 if (!(data->bound && data->bound[idof] >= DIRICHLET)) {
1284 jdof = data->local_dof[i];
1285 for (k = 0; k < N_VERTICES(dim); k++)
1286 AXPY_DOW(data->ipol[jdof][k], h[data->dof_parent[idof][k]], h[idof]);
1287 }
1288 }
1289 }
1290 return;
1291 }
1292
1293 /****************************************************************************/
1294
get_BPX_precon_s(const DOF_MATRIX * matrix,const DOF_SCHAR_VEC * bound,int info)1295 static const PRECON *get_BPX_precon_s(const DOF_MATRIX *matrix,
1296 const DOF_SCHAR_VEC *bound,
1297 int info)
1298 {
1299 PRECON *precon;
1300 HB_DATA *data;
1301 SCRATCH_MEM scratch;
1302 const FE_SPACE *fe_space = matrix->row_fe_space;
1303
1304 if (bound && !FE_SPACE_EQ_P(bound->fe_space, fe_space)) {
1305 ERROR("different fe spaces ?\n");
1306 return(NULL);
1307 }
1308
1309 #if 0
1310 /* this stuff does not work with DOF-space holes */
1311 dof_compress(fe_space->mesh);
1312 #endif
1313
1314 SCRATCH_MEM_INIT(scratch);
1315 data = SCRATCH_MEM_CALLOC(scratch, 1, HB_DATA);
1316 SCRATCH_MEM_CPY(data->scratch, scratch);
1317
1318 precon = &data->precon;
1319
1320 data->matrix = matrix;
1321 data->fe_space = fe_space;
1322 data->bound_dv = bound;
1323 BNDRY_FLAGS_CPY(data->dirichlet_bndry, matrix->dirichlet_bndry);
1324
1325 data->info = info;
1326
1327 precon->precon_data = data;
1328 precon->init_precon = init_BPX_precon;
1329 precon->precon = BPX_precon_s;
1330 precon->exit_precon = exit_BPX_precon;
1331
1332 return(precon);
1333 }
1334
get_BPX_precon_d(const DOF_MATRIX * matrix,const DOF_SCHAR_VEC * bound,int info)1335 static const PRECON *get_BPX_precon_d(const DOF_MATRIX *matrix,
1336 const DOF_SCHAR_VEC *bound,
1337 int info)
1338 {
1339 PRECON *precon;
1340 HB_DATA *data;
1341 SCRATCH_MEM scratch;
1342 const FE_SPACE *fe_space = matrix->row_fe_space;
1343
1344 if (bound && !FE_SPACE_EQ_P(bound->fe_space, fe_space)) {
1345 ERROR("different fe spaces ?\n");
1346 return(NULL);
1347 }
1348
1349 #if 0
1350 /* this stuff does not work with DOF-space holes */
1351 dof_compress(fe_space->mesh);
1352 #endif
1353
1354 SCRATCH_MEM_INIT(scratch);
1355 data = SCRATCH_MEM_CALLOC(scratch, 1, HB_DATA);
1356 SCRATCH_MEM_CPY(data->scratch, scratch);
1357
1358 precon = &data->precon;
1359
1360 data->matrix = matrix;
1361 data->fe_space = fe_space;
1362 data->bound_dv = bound;
1363 BNDRY_FLAGS_CPY(data->dirichlet_bndry, matrix->dirichlet_bndry);
1364
1365 data->info = info;
1366
1367 precon->precon_data = data;
1368 precon->init_precon = init_BPX_precon;
1369 precon->precon = BPX_precon_d;
1370 precon->exit_precon = exit_BPX_precon;
1371
1372 return precon;
1373 }
1374
1375 /****************************************************************************/
1376
get_HB_precon(const DOF_MATRIX * matrix,const DOF_SCHAR_VEC * bound,int info)1377 const PRECON *get_HB_precon(const DOF_MATRIX *matrix,
1378 const DOF_SCHAR_VEC *bound,
1379 int info)
1380 {
1381 const FE_SPACE *fe_space = matrix->row_fe_space;
1382
1383 TEST_EXIT(fe_space->bas_fcts->rdim == 1,
1384 "This cannot work for exotic FE-spaces.\n");
1385 if (fe_space->rdim == 1) {
1386 return get_HB_precon_s(matrix, bound, info);
1387 } else {
1388 return get_HB_precon_d(matrix, bound, info);
1389 }
1390 }
1391
get_BPX_precon(const DOF_MATRIX * matrix,const DOF_SCHAR_VEC * bound,int info)1392 const PRECON *get_BPX_precon(const DOF_MATRIX *matrix,
1393 const DOF_SCHAR_VEC *bound,
1394 int info)
1395 {
1396 const FE_SPACE *fe_space = matrix->row_fe_space;
1397
1398 TEST_EXIT(fe_space->bas_fcts->rdim == 1,
1399 "This cannot work for exotic FE-spaces.\n");
1400 if (fe_space->rdim == 1) {
1401 return get_BPX_precon_s(matrix, bound, info);
1402 } else {
1403 return get_BPX_precon_d(matrix, bound, info);
1404 }
1405 }
1406
1407
1408