1 /*******************************************************************************
2 *
3 * ALBERTA: an Adaptive multi Level finite element toolbox using
4 * Bisectioning refinement and Error control by Residual
5 * Techniques
6 *
7 * www.alberta-fem.de
8 *
9 ******************************************************************************
10 *
11 * File: quasi-stokes.c
12 *
13 * Description: solver for a simple Quasi-Stokes problem:
14 *
15 * \mu u - \nu\Delta u + \nabla p = f in \Omega
16 * \nabla\cdot u = 0 in \Omega
17 * \sigma(u, p) n = F on \partial \Omega
18 *
19 * with \sigma(u, p) = \nu D(u) - I p, D(u) = \nabla u + (\nable u)^t, n
20 * is the outer normal to \partial\Omega.
21 *
22 * To incorporate the stress boundary-conditions it is necessary to
23 * use the symmetric gradient in the distributional formulation, the
24 * resulting stiffness-matrix is a DIM_OF_WORLD x DIM_OF_WORLD block
25 * matrix. The physical meaning of the boundary conditions is that the
26 * intrinsic stresses of the fluid are counter-balanced by the
27 * force-density F on the boundary.
28 *
29 * The code below can actually also be used to solve a Quasi-Stokes
30 * problem with Dirichlet boundary conditions; there is a
31 * corresponding parameter in the init-file.
32 *
33 ******************************************************************************
34 *
35 * author(s): Claus-Justus Heine
36 * Abteilung fuer Angewandte Mathematik
37 * Albert-Ludwigs-Universitaet Freiburg
38 * Hermann-Herder-Str. 10
39 * 79104 Freiburg
40 * Germany
41 * Claus.Heine@Mathematik.Uni-Freiburg.DE
42 *
43 * (c) by C.-J. Heine (2007-2012)
44 *
45 ******************************************************************************/
46
47 #include <alberta/alberta.h>
48 #include <alberta/albas.h>
49
50 #include "alberta-demo.h"
51
52 static int do_graphics = true;
53
54 /*******************************************************************************
55 * global variables: finite element space, discrete solution
56 * load vector and system matrix
57 *
58 * These variables are kept global because they are shared across build(),
59 * solve() and estimate().
60 ******************************************************************************/
61
62 static const FE_SPACE *u_fe_space; /* velocity FE-space */
63 static const FE_SPACE *p_fe_space; /* pressure FE-space */
64 static DOF_REAL_VEC_D *u_h; /* discrete velocity */
65 static DOF_REAL_VEC *p_h; /* discrete pressure */
66
67 static DOF_REAL_VEC_D *f_h; /* RHS for velocity */
68 static DOF_REAL_VEC *g_h; /* RHS for pressure (if necessary) */
69 static DOF_SCHAR_VEC *bound; /* boundary vector */
70 static DOF_MATRIX *A; /* system matrix */
71
72 static DOF_MATRIX *B; /* discrete pressure gradient */
73 static DOF_MATRIX *Bt; /* discrete velocity divergence */
74 static DOF_MATRIX *Yproj; /* projection/precon matrix for B^\ast */
75 static DOF_MATRIX *Yprec; /* projection/precon matrix for B^\ast */
76
77 /* Quasi-Stokes parameters */
78 static REAL mu = 1.0;
79 static REAL nu = 1.0;
80 static REAL alpha_robin = 1.0;
81
82 /* Parameters for the exact "solution" */
83 static REAL gauss_bell = 10.0; /* factor for the pressure */
84 static REAL_D p_shift; /* translation for the pressure */
85 static REAL_D v_shift; /* translation for the velocity */
86
87 static BNDRY_FLAGS dirichlet_bndry; /* Dirichlet part of the bndry, if any. */
88 static BNDRY_FLAGS stress_bndry; /* Complement of dirichlet_mask. */
89
90 /*******************************************************************************
91 * struct ellipt_leaf_data: structure for storing one REAL value on each
92 * leaf element as LEAF_DATA
93 * rw_el_est(): return a pointer to the memory for storing the element
94 * estimate (stored as LEAF_DATA), called by ellipt_est()
95 * get_el_est(): return the value of the element estimates (from LEAF_DATA),
96 * called by adapt_method_stat() and graphics()
97 ******************************************************************************/
98
99 struct ellipt_leaf_data
100 {
101 REAL estimate; /* one real for the estimate */
102 };
103
rw_el_est(EL * el)104 static REAL *rw_el_est(EL *el)
105 {
106 if (IS_LEAF_EL(el))
107 return &((struct ellipt_leaf_data *)LEAF_DATA(el))->estimate;
108 else
109 return NULL;
110 }
111
get_el_est(EL * el)112 static REAL get_el_est(EL *el)
113 {
114 if (IS_LEAF_EL(el))
115 return ((struct ellipt_leaf_data *)LEAF_DATA(el))->estimate;
116 else
117 return 0.0;
118 }
119
120 /*******************************************************************************
121 * For test purposes: exact solution and its gradient (optional)
122 ******************************************************************************/
123
124 /*
125 * u(x) = \frac{x}{|x|^d} (make sure 0 lives outside of the domain ...
126 */
u(const REAL_D _x,REAL_D res)127 static const REAL *u(const REAL_D _x, REAL_D res)
128 {
129 static REAL_D space;
130 REAL r;
131
132 if (!res) {
133 res = space;
134 }
135
136 AXPBY_DOW(1.0, _x, -1.0, v_shift, res);
137
138 r = NORM_DOW(res);
139
140 SCAL_DOW(1.0/POW_DOW(r), res);
141
142 return res;
143 }
144
145 /*
146 * [(r^2-x^2*d)*r^(-d-2), -d*x*y*r^(-d-2), -d*x*z*r^(-d-2)]
147 * etc.
148 */
grd_u(const REAL_D _x,REAL_DD res)149 static const REAL_D *grd_u(const REAL_D _x, REAL_DD res)
150 {
151 static REAL_DD space;
152 REAL_D x;
153 REAL r;
154 int i;
155
156 if (!res) {
157 res = space;
158 }
159
160 AXPBY_DOW(1.0, _x, -1.0, v_shift, x);
161
162 r = NORM_DOW(x);
163
164 for (i = 0; i < DIM_OF_WORLD; i++) {
165 AXEY_DOW(-(REAL)DIM_OF_WORLD*x[i], x, res[i]);
166 res[i][i] += SQR(r);
167 SCAL_DOW(pow(r, (REAL)(-DIM_OF_WORLD - 2)), res[i]);
168 }
169 return (const REAL_D *)res;
170 }
171
p(const REAL_D _x)172 static REAL p(const REAL_D _x)
173 {
174 REAL_D x;
175
176 AXPBY_DOW(1.0, _x, -1.0, p_shift, x);
177
178 return exp(-gauss_bell*SCP_DOW(x,x));
179 }
180
grd_p(const REAL_D _x,REAL_D res)181 static const REAL *grd_p(const REAL_D _x, REAL_D res)
182 {
183 static REAL_D space;
184 REAL_D x;
185
186 if (!res) {
187 res = space;
188 }
189
190 AXPBY_DOW(1.0, _x, -1.0, p_shift, x);
191
192 AXEY_DOW(-2.0*gauss_bell*p(_x), x, res);
193
194 return res;
195 }
196
197 /*******************************************************************************
198 * problem data: right hand side, boundary values
199 ******************************************************************************/
200
201 /* Dirichlet boundary values */
g(const REAL_D x,REAL_D res)202 static const REAL *g(const REAL_D x, REAL_D res)
203 {
204 return u(x, res);
205 }
206
207 /* \mu u - \nu \Delta u + \nabla p */
f(const REAL_D _x,REAL_D res)208 static const REAL *f(const REAL_D _x, REAL_D res)
209 {
210 static REAL_D space;
211 REAL_D u_val;
212
213 if (!res) {
214 res = space;
215 }
216
217 /* u is harmonic */
218 AXPY_DOW(mu, u(_x, u_val), (REAL *)grd_p(_x, res));
219
220 return res;
221 }
222
223 /* stress tensor on the boundary, needed for non-zero stress boundary
224 * conditions.
225 */
g_stress(const REAL_D x,const REAL_D normal,REAL_D res)226 static const REAL *g_stress(const REAL_D x,
227 const REAL_D normal,
228 REAL_D res)
229 {
230 static REAL_D space;
231 int i;
232 REAL_DD Dv;
233 REAL p_val;
234
235 if (!res) {
236 res = space;
237 }
238
239 grd_u(x, Dv); /* in our example grd_u is already symmetric */
240 p_val = p(x);
241 MSCAL_DOW(2.0*nu, Dv);
242 for (i = 0; i < DIM_OF_WORLD; i++) {
243 Dv[i][i] -= p_val;
244 }
245
246 SET_DOW(0.0, res);
247 MV_DOW((const REAL_D *)Dv, normal, res);
248
249 return res;
250 }
251
252 /*******************************************************************************
253 * build(): assemblage of the linear system: matrix, load vector,
254 * boundary values, called by adapt_method_stat()
255 * on the first call initialize u_h, f_h, matrix and information
256 * for assembling the system matrix
257 *
258 * struct op_info: structure for passing information from init_element() to
259 * LALt()
260 *
261 * u_init_element(), p_init_element():
262 * initialization on the element; calculates the
263 * coordinates and |det DF_S| used by LALt; passes these
264 * values to LALt via user_data,
265 * called on each element by update_matrix()
266 *
267 * LDDLt(): implementation of Lambda DD Lambda^t for the -Delta v
268 * term in the stress-tensor formulation
269 *
270 * LALt(): implementation of Lambda id Lambda^t for -Delta p,
271 * called by update_matrix() after init_element()
272 *
273 * B_Lb1(): operator kernel (just Lambda) for \phi\nabla\psi
274 *
275 * c(): mass matrix for the projection to the pressure space
276 *
277 ******************************************************************************/
278
279 typedef struct op_data
280 {
281 REAL_BD Lambda; /* the gradient of the barycentric coordinates */
282 REAL det; /* |det D F_S| */
283 REAL LALt_factor; /* a factor for the 2nd order element matrices */
284 REAL c_factor; /* a factor for the 0th order element matrices */
285 REAL wall_det[N_WALLS_MAX];
286 } OP_DATA;
287
u_init_element(const EL_INFO * el_info,const QUAD * quad[3],void * ud)288 static bool u_init_element(const EL_INFO *el_info,
289 const QUAD *quad[3],
290 void *ud)
291 {
292 OP_DATA *data = (OP_DATA *)ud;
293
294 data->det = el_grd_lambda(el_info, data->Lambda);
295
296 data->LALt_factor = data->det * nu;
297 data->c_factor = data->det * mu;
298
299 /* Return false, i.e. not parametric */
300 return false;
301 }
302
p_init_element(const EL_INFO * el_info,const QUAD * quad[3],void * ud)303 static bool p_init_element(const EL_INFO *el_info,
304 const QUAD *quad[3],
305 void *ud)
306 {
307 OP_DATA *data = (OP_DATA *)ud;
308
309 data->det = el_grd_lambda(el_info, data->Lambda);
310
311 data->LALt_factor = data->det;
312 data->c_factor = data->det;
313
314 return 0;
315 }
316
dg_wall_init_element(const EL_INFO * el_info,int wall,const WALL_QUAD * quad[3],void * ud)317 static bool dg_wall_init_element(const EL_INFO *el_info, int wall,
318 const WALL_QUAD *quad[3], void *ud)
319 {
320 OP_DATA *info = (OP_DATA *)ud;
321
322 info->wall_det[wall] = 1e3;
323
324 info->c_factor = info->wall_det[wall];
325
326 return false;
327 }
328
fv_wall_init_element(const EL_INFO * el_info,int wall,const WALL_QUAD * quad[3],void * ud)329 static bool fv_wall_init_element(const EL_INFO *el_info, int wall,
330 const WALL_QUAD *quad[3], void *ud)
331 {
332 OP_DATA *info = (OP_DATA *)ud;
333 REAL det = get_wall_normal(el_info, wall, NULL);
334 /* Use a very simplistic finite-volume method */
335 #define VAL (1.0/(REAL)N_LAMBDA_MAX)
336 const REAL_B c_b = INIT_BARY_MAX(VAL, VAL, VAL, VAL);
337 #undef VAL
338 EL_INFO neigh_info[1];
339 const EL_GEOM_CACHE *elgc =
340 fill_el_geom_cache(el_info, FILL_EL_WALL_REL_ORIENTATION(wall));
341 REAL_D c, cn;
342 REAL delta;
343
344 fill_neigh_el_info(neigh_info, el_info, wall, elgc->rel_orientation[wall]);
345
346 coord_to_world(el_info, c_b, c);
347 coord_to_world(neigh_info, c_b, cn);
348 delta = DIST_DOW(c, cn);
349
350 info->c_factor = info->wall_det[wall] = det / delta;
351
352 return false;
353 }
354
dg_fv_neigh_init_element(const EL_INFO * el_info,int wall,const WALL_QUAD * quad[3],void * ud)355 static bool dg_fv_neigh_init_element(const EL_INFO *el_info, int wall,
356 const WALL_QUAD *quad[3], void *ud)
357 {
358 OP_DATA *info = (OP_DATA *)ud;
359
360 info->c_factor = -info->wall_det[wall];
361
362 return false;
363 }
364
LALt(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)365 static const REAL_B *LALt(const EL_INFO *el_info,
366 const QUAD *quad, int iq, void *ud)
367 {
368 OP_DATA *data = (OP_DATA *)ud;
369 int i, j, dim = el_info->mesh->dim;
370 static REAL_BB LALt;
371
372 for (i = 0; i <= dim; i++) {
373 for (j = i; j <= dim; j++) {
374 LALt[j][i] = LALt[i][j] =
375 data->LALt_factor*SCP_DOW(data->Lambda[i], data->Lambda[j]);
376 }
377 }
378
379 return (const REAL_B *)LALt;
380 }
381
382 /* Deformation tensor. This is 1/2 D(\phi^l_k):D(\phi^n_m) transformed
383 * to the reference element, where the components of the vector-valued
384 * test functions are given by (\phi^l_k)_j = \delta_{kj}\phi^l for an
385 * "ordinary" scalar basis function.
386 *
387 * We have to assemble the block-matrices
388 * ${{\cal B}^{\alpha\beta}}\in\R^{d\times d}$:
389 *
390 * ({{\cal B}^{\alpha\beta}})_{km}
391 * =
392 * \delta_{km}\,\Lambda^{\alpha r}\,\Lambda^{\beta r}
393 * +
394 * \Lambda^{\alpha m}\,\Lambda^{\beta k}.
395 */
396 #if 0
397 static const REAL_BDD *LDDLt(const EL_INFO *el_info, const QUAD *quad,
398 int iq, void *ud)
399 {
400 static REAL_BBDD LDDLt;
401 OP_DATA *info = (OP_DATA *)ud;
402 int m, n, i, j, dim = el_info->mesh->dim;
403 REAL factor = info->LALt_factor;
404
405 for (m = 0; m < DIM_OF_WORLD; m++)
406 for (n = 0; n < DIM_OF_WORLD; n++)
407 for (i = 0; i <= dim; i++)
408 for (j = 0; j <= dim; j++) {
409 LDDLt[i][j][m][n] = 0.0;
410
411 if(m == n) {
412 LDDLt[i][j][m][m] += SCP_DOW(info->Lambda[i], info->Lambda[j]);
413 }
414
415 LDDLt[i][j][m][n] += info->Lambda[i][n] * info->Lambda[j][m];
416 LDDLt[i][j][m][n] *= factor;
417 }
418 return (const REAL_BDD *)LDDLt;
419 }
420 #else
LDDLt(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)421 static const REAL_BDD *LDDLt(const EL_INFO *el_info,
422 const QUAD *quad, int iq, void *ud)
423 {
424 static REAL_BBDD LDDLt;
425 OP_DATA *info = (OP_DATA *)ud;
426 REAL value;
427 int alpha, beta, k, m;
428 int dim = el_info->mesh->dim;
429 REAL factor = info->LALt_factor;
430
431 for (alpha = 0; alpha < N_VERTICES(dim); alpha++) {
432 /* The first part is just LLt reordered such that it matches the
433 * block-matrix conventions
434 */
435 MSET_DOW(0.0, LDDLt[alpha][alpha]);
436 value = factor*SCP_DOW(info->Lambda[alpha], info->Lambda[alpha]);
437 MSCMAXPY_DOW(1.0, value, LDDLt[alpha][alpha]);
438 for (beta = alpha+1; beta < N_VERTICES(dim); beta++) {
439 MSET_DOW(0.0, LDDLt[alpha][beta]);
440 MSET_DOW(0.0, LDDLt[beta][alpha]);
441 value = factor*SCP_DOW(info->Lambda[alpha], info->Lambda[beta]);
442 MSCMAXPY_DOW(1.0, value, LDDLt[alpha][beta]);
443 MSCMAXPY_DOW(1.0, value, LDDLt[beta][alpha]);
444 }
445
446 /* The second part is symmtric w.r.t. to the exchange of index
447 * pairs. We do loop unrolling here to reduce the number of loads
448 * required to store stuff into the 4 dimensional array.
449 */
450 for (m = 0; m < DIM_OF_WORLD; m++) {
451 /* (alpha, m, alpha, m) */
452 value = factor*SQR(info->Lambda[alpha][m]);
453 LDDLt[alpha][alpha][m][m] += value;
454 for (k = m+1; k < DIM_OF_WORLD; k++) {
455 /* (alpha, m, alpha, k)
456 * (alpha, k, alpha, m)
457 */
458 value = factor*info->Lambda[alpha][m]*info->Lambda[alpha][k];
459 LDDLt[alpha][alpha][m][k] += value;
460 LDDLt[alpha][alpha][k][m] += value;
461 }
462 for (beta = alpha+1; beta < N_VERTICES(dim); beta++) {
463 /* (alpha, m, beta, m)
464 * (beta, m, alpha, m)
465 */
466 value = factor*info->Lambda[alpha][m]*info->Lambda[beta][m];
467 LDDLt[alpha][beta][m][m] += value;
468 LDDLt[beta][alpha][m][m] += value;
469 for (k = m+1; k < DIM_OF_WORLD; k++) {
470 /* (alpha, m, beta, k)
471 * (beta, k, alpha, m)
472 */
473 value = factor*info->Lambda[alpha][m]*info->Lambda[beta][k];
474 LDDLt[alpha][beta][k][m] += value;
475 LDDLt[beta][alpha][m][k] += value;
476 value = factor*info->Lambda[alpha][k]*info->Lambda[beta][m];
477 LDDLt[alpha][beta][m][k] += value;
478 LDDLt[beta][alpha][k][m] += value;
479 }
480 }
481 }
482 }
483 return (const REAL_BDD *)LDDLt;
484 }
485 #endif
486
487 /* The first order term for the divergence constraint */
B_Lb1(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)488 const REAL_D *B_Lb1(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
489 {
490 OP_DATA *info = (OP_DATA *)ud;
491 static REAL_BD B_Lb1;
492 int i;
493
494 for (i = 0; i < N_VERTICES(el_info->mesh->dim); i++) {
495 AXEY_DOW(-info->det, info->Lambda[i], B_Lb1[i]);
496 }
497
498 return (const REAL_D *)B_Lb1;
499 }
500
c(const EL_INFO * el_info,const QUAD * quad,int iq,void * ud)501 static REAL c(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
502 {
503 OP_DATA *info = (OP_DATA *)ud;
504
505 return info->c_factor;
506 }
507
build(MESH * mesh,U_CHAR flag)508 static void build(MESH *mesh, U_CHAR flag)
509 {
510 FUNCNAME("build");
511 static const EL_MATRIX_INFO *matrix_info;
512 static const EL_MATRIX_INFO *Yproj_minfo;
513 static const EL_MATRIX_INFO *Yprec_minfo;
514 static EL_MATRIX_INFO B_minfo;
515 static OP_DATA opi_udata;
516
517 const QUAD *quad;
518 REAL flux;
519
520 dof_compress(mesh);
521 MSG("%d DOFs for %s\n", dof_real_vec_d_length(u_fe_space), u_fe_space->name);
522 MSG("%d DOFs for %s\n", p_fe_space->admin->size_used, p_fe_space->name);
523
524 if (!matrix_info) { /* information for matrix assembling */
525 OPERATOR_INFO o_info = { NULL, };
526 OPERATOR_INFO B_oinfo = { NULL, };
527 OPERATOR_INFO Yproj_oinfo = { NULL, };
528 OPERATOR_INFO Yprec_oinfo = { NULL, };
529 bool p_disc;
530
531 o_info.row_fe_space = o_info.col_fe_space = u_fe_space;
532 o_info.init_element = u_init_element;
533 o_info.user_data = &opi_udata;
534 o_info.LALt.real_dd = LDDLt;
535 o_info.LALt_type = MATENT_REAL_DD;
536 o_info.LALt_pw_const = true; /* pw const. assemblage is faster */
537 o_info.LALt_symmetric = false /*true*/; /* symmetric assemblage is faster */
538 o_info.c.real = c;
539 o_info.c_type = MATENT_REAL;
540 o_info.c_pw_const = true;
541 BNDRY_FLAGS_CPY(o_info.dirichlet_bndry,
542 dirichlet_bndry); /* Dirichlet bc on part of the bndry? */
543 o_info.fill_flag = CALL_LEAF_EL|FILL_COORDS;
544 matrix_info = fill_matrix_info(&o_info, NULL);
545
546 /* Operator-info for divergence constraint, only the B operator,
547 * B^\ast will be generated by transposition of B and projection
548 * to the pressure space.
549 */
550 B_oinfo.row_fe_space = u_fe_space;
551 B_oinfo.col_fe_space = p_fe_space;
552 B_oinfo.init_element = p_init_element;
553 B_oinfo.user_data = o_info.user_data;
554 B_oinfo.Lb1.real_d = B_Lb1;
555 B_oinfo.Lb_type = MATENT_REAL_D;
556 B_oinfo.Lb1_pw_const = true;
557 B_oinfo.fill_flag = CALL_LEAF_EL|FILL_COORDS;
558 fill_matrix_info(&B_oinfo, &B_minfo);
559
560 /* Projection to the pressure space: B^\ast = Pr_Y \circ B^{tr} */
561 Yproj_oinfo.row_fe_space = Yproj_oinfo.col_fe_space = p_fe_space;
562 Yproj_oinfo.init_element = p_init_element;
563 Yproj_oinfo.user_data = o_info.user_data;
564 Yproj_oinfo.c.real = c;
565 Yproj_oinfo.c_pw_const = true;
566 Yproj_oinfo.fill_flag = CALL_LEAF_EL|FILL_COORDS;
567 Yproj_minfo = fill_matrix_info(&Yproj_oinfo, NULL);
568
569 /* (Part of) the preconditioner for the Schur-Complement operator;
570 * the Poissoin-part of the pre-conditioner corresponds to the
571 * mass-part of the Quasi-Stokes operator.
572 */
573 Yprec_oinfo.row_fe_space = Yprec_oinfo.col_fe_space = p_fe_space;
574 Yprec_oinfo.init_element = p_init_element;
575 Yprec_oinfo.user_data = o_info.user_data;
576 Yprec_oinfo.LALt.real = LALt;
577 Yprec_oinfo.LALt_symmetric = true;
578 Yprec_oinfo.LALt_pw_const = true;
579 Yprec_oinfo.fill_flag = CALL_LEAF_EL|FILL_COORDS;
580
581 p_disc =
582 p_fe_space->bas_fcts->n_dof[CENTER] == p_fe_space->bas_fcts->n_bas_fcts;
583 if (!p_disc) {
584 Yprec_minfo = fill_matrix_info(&Yprec_oinfo, NULL);
585 } else {
586 /* Discontinuous pressure, use a finite volume (degree == 0) or
587 * DG method (degree > 0) to solve the Poisson problem.
588 */
589 BNDRY_OPERATOR_INFO el_bop_info = { NULL, };
590 BNDRY_OPERATOR_INFO neigh_bop_info = { NULL, };
591 bool do_fv = p_fe_space->bas_fcts->degree == 0;
592
593 /* ALBERTA requires that the boundary contributions are split
594 * into one part which is assembled using only the current's
595 * element local basis functions, and a second part which pairs
596 * the given element with its neighbour. In our case both
597 * contributions are simple mass-matrices
598 */
599 el_bop_info.row_fe_space = p_fe_space;
600 el_bop_info.init_element =
601 do_fv ? fv_wall_init_element : dg_wall_init_element;
602 el_bop_info.c.real = c;
603 el_bop_info.c_pw_const = true;
604 el_bop_info.user_data = o_info.user_data;
605 BNDRY_FLAGS_INIT(el_bop_info.bndry_type); /* all interior edges */
606
607 el_bop_info.fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_OPP_COORDS;
608
609 /* The contribution of the pairing with the neighbouring basis
610 * functions is almost the same:
611 */
612 neigh_bop_info = el_bop_info;
613 neigh_bop_info.init_element = dg_fv_neigh_init_element;
614 neigh_bop_info.discontinuous = true;
615
616 Yprec_minfo = fill_matrix_info_ext(
617 NULL, do_fv ? NULL : &Yprec_oinfo, &el_bop_info, &neigh_bop_info, NULL);
618 }
619 }
620
621 /* assembling of matrix */
622 clear_dof_matrix(A);
623 update_matrix(A, matrix_info, NoTranspose);
624
625 /* assembling of load vector */
626 dof_set_dow(0.0, f_h);
627 quad = get_quadrature(mesh->dim, 2*u_fe_space->bas_fcts->degree - 2);
628 L2scp_fct_bas_dow(f, quad, f_h);
629
630 /* boundary values, we have to fill in the "bound" vector for
631 * sp_dirichlet_bound_sd() below, and add a boundary integral for
632 * inhomogeneous stress boundary conditions.
633 *
634 * For a pure Stokes problem (mu == 0.0) we also request an
635 * automatic mean-value correction of the right-hand side; this is
636 * triggered by "bndry_type == NEUMANN && robin_alpha < 0.0".
637 */
638 boundary_conditions_dow(NULL /* matrix */, f_h, u_h, bound,
639 dirichlet_bndry,
640 g, g_stress,
641 mu == 0.0 ? -1.0 : 0.0 /* robin_alpha */,
642 NULL);
643
644 /* Now go for the divergence constraint */
645
646 /* First for B */
647 clear_dof_matrix(B);
648 BNDRY_FLAGS_CPY(B_minfo.dirichlet_bndry, dirichlet_bndry);
649 update_matrix(B, &B_minfo, NoTranspose);
650
651 /* Then for Bt. Ideally, one would loop just once over the mesh
652 * and assemble both matrices together.
653 */
654 clear_dof_matrix(Bt);
655 BNDRY_FLAGS_INIT(B_minfo.dirichlet_bndry);
656 update_matrix(Bt, &B_minfo, Transpose);
657
658 /* In order to have compatible boundary conditions on the discrete
659 * level we have to make sure that either the boundary integral over
660 * the boundary values vanishes (discretely), or we have to perturb
661 * the RHS for the pressure, i.e. solve a (slightly) inhomogeneous
662 * saddle-point problem. This is what we are doing here because it
663 * is conceptual simpler.
664 */
665 dof_set(0.0, g_h);
666 flux = sp_dirichlet_bound_dow_scl(NoTranspose, Bt, bound, u_h, g_h);
667
668 MSG("Total flux-excess: %e\n", flux);
669
670 /* Finally the projection/preconditioner for the pressure space */
671
672 /* Assemble the projection; at the same time the part of the
673 * preconditioner "responsible" for the 2nd-order part of the
674 * Quasi-Stokes operator.
675 */
676 clear_dof_matrix(Yproj);
677 update_matrix(Yproj, Yproj_minfo, NoTranspose);
678
679 /* Assemble the part of the preconditiner "responsible" for the
680 * 0th-order part of the Quasi-Stokes operator. We need Robin
681 * boundary conditions because of the Neumann-style boundary
682 * conditions in the velocity space.
683 */
684 clear_dof_matrix(Yprec);
685 update_matrix(Yprec, Yprec_minfo, NoTranspose);
686 robin_bound(Yprec, stress_bndry, alpha_robin, NULL /* bndry_quad */, 1.0);
687
688 /* That's it */
689 }
690
691 /*******************************************************************************
692 * solve(): solve the linear system, called by adapt_method_stat()
693 ******************************************************************************/
694
solve(MESH * mesh)695 static void solve(MESH *mesh)
696 {
697 FUNCNAME("solve");
698 static REAL tol = 1.e-8, tol_incr = 1e2;
699 static int miter = 1000, info = 2, restart = 0;
700 static OEM_SOLVER solver = NoSolver, A_solver = NoSolver;
701 static OEM_SOLVER Yproj_solver = NoSolver, Yprec_solver = NoSolver;
702 static int A_miter, A_restart;
703 static int Yproj_miter, Yprec_miter;
704 const PRECON *A_prec, *Yproj_prec, *Yprec_prec;
705 static PRECON_TYPE A_prec_type = {
706 BlkDiagPrecon,
707 };
708 static PRECON_TYPE Yproj_prec_type = {
709 DiagPrecon,
710 };
711 static PRECON_TYPE Yprec_prec_type = {
712 DiagPrecon,
713 };
714
715 if (solver == NoSolver) {
716 struct __precon_type *first = &A_prec_type.param.BlkDiag.precon[0];
717 int i;
718
719 GET_PARAMETER(1, "sp->solver", "%d", &solver);
720 GET_PARAMETER(1, "sp->solver tolerance", "%f", &tol);
721 GET_PARAMETER(1, "sp->solver tol incr", "%f", &tol_incr);
722 GET_PARAMETER(1, "sp->solver max iteration", "%d", &miter);
723 GET_PARAMETER(1, "sp->solver info", "%d", &info);
724 if (solver == GMRes) {
725 GET_PARAMETER(1, "sp->solver restart", "%d", &restart);
726 }
727
728 GET_PARAMETER(1, "sp->Auf solver", "%d", &A_solver);
729 GET_PARAMETER(1, "sp->Auf max iteration", "%d", &A_miter);
730 GET_PARAMETER(1, "sp->Auf precon", "%d", &first->type);
731 for (i = 1; i < N_BLOCK_PRECON_MAX; ++i)
732 A_prec_type.param.BlkDiag.precon[i].type = DiagPrecon;
733
734 if (first->type == __SSORPrecon) {
735 GET_PARAMETER(1, "sp->Auf precon ssor omega", "%f",
736 &first->param.__SSOR.omega);
737 GET_PARAMETER(1, "sp->Auf precon ssor iter", "%d",
738 &first->param.__SSOR.n_iter);
739 }
740 if (first->type == ILUkPrecon)
741 GET_PARAMETER(1, "sp->Auf precon ilu(k)", "%d",
742 &first->param.ILUk.level);
743
744
745 if (A_solver == GMRes) {
746 GET_PARAMETER(1, "sp->Auf restart", "%d", &A_restart);
747 }
748
749 GET_PARAMETER(1, "sp->Yproj solver", "%d", &Yproj_solver);
750 GET_PARAMETER(1, "sp->Yproj max iteration", "%d", &Yproj_miter);
751 GET_PARAMETER(1, "sp->Yproj precon", "%d", &Yproj_prec_type.type);
752 if (Yproj_prec_type.type == __SSORPrecon) {
753 GET_PARAMETER(1, "sp->Yproj precon ssor omega", "%f",
754 &Yproj_prec_type.param.__SSOR.omega);
755 GET_PARAMETER(1, "sp->Yproj precon ssor iter", "%d",
756 &Yproj_prec_type.param.__SSOR.n_iter);
757 }
758 if (Yproj_prec_type.type == ILUkPrecon)
759 GET_PARAMETER(1, "sp->Yproj precon ilu(k)", "%d",
760 &Yproj_prec_type.param.ILUk.level);
761
762
763
764 GET_PARAMETER(1, "sp->Yprec solver", "%d", &Yprec_solver);
765 GET_PARAMETER(1, "sp->Yprec max iteration", "%d", &Yprec_miter);
766 GET_PARAMETER(1, "sp->Yprec precon", "%d", &Yprec_prec_type.type);
767 if (Yprec_prec_type.type == __SSORPrecon) {
768 GET_PARAMETER(1, "sp->Yprec precon ssor omega", "%f",
769 &Yprec_prec_type.param.__SSOR.omega);
770 GET_PARAMETER(1, "sp->Yprec precon ssor iter", "%d",
771 &Yprec_prec_type.param.__SSOR.n_iter);
772 }
773 if (Yprec_prec_type.type == ILUkPrecon)
774 GET_PARAMETER(1, "sp->Yprec precon ilu(k)", "%d",
775 &Yprec_prec_type.param.ILUk.level);
776 }
777
778 A_prec = init_precon_from_type(A, NULL, MAX(0, info-3), &A_prec_type);
779 Yproj_prec = init_precon_from_type(Yproj, NULL /* bound */, MAX(0, info-3),
780 &Yproj_prec_type);
781 Yprec_prec = init_precon_from_type(Yprec, NULL /* bound */, MAX(0, info-3),
782 &Yprec_prec_type);
783
784 oem_sp_solve_dow_scl(solver, tol, tol_incr, miter, info,
785 A, NULL /* bound */, A_solver, A_miter, A_prec,
786 B, Bt,
787 Yproj, Yproj_solver, Yproj_miter, Yproj_prec,
788 mu > 0 ? Yprec : NULL, Yprec_solver, Yprec_miter,
789 Yprec_prec,
790 nu /* Yproj_frac */, mu /* Yprec_frac */,
791 f_h, g_h, u_h, p_h);
792
793 }
794
795 /*******************************************************************************
796 * Functions for error estimate:
797 *
798 * estimate(): compute a residual error estimate and the exact error
799 * (the latter for testing purposes). The estimate is performed for
800 * the velocity only, the pressure enters as lower-order term on the
801 * RHS.
802 *
803 * r(): calculates the lower order terms of the element residual on
804 * each element at the quadrature node iq of quad argument to
805 * ellipt_est() and called by ellipt_est()
806 *
807 * stress_res(): residual of the stress boundary conditions.
808 ******************************************************************************/
809
810 /* Lower order terms: right-hand side, the 0th order term, the
811 * gradient of the discrete pressure.
812 */
r(REAL_D fx,const EL_INFO * el_info,const QUAD * quad,int iq,const REAL_D uh_at_qp,const REAL_DD grd_uh_at_qp)813 static const REAL *r(REAL_D fx,
814 const EL_INFO *el_info,
815 const QUAD *quad,
816 int iq,
817 const REAL_D uh_at_qp, const REAL_DD grd_uh_at_qp)
818 {
819 static REAL_D space;
820 static const QUAD_FAST *qfast;
821 static EL *el;
822 static EL_REAL_VEC *ph_loc;
823 static const EL_GEOM_CACHE *elgc;
824 static const QUAD_EL_CACHE *qelc;
825 REAL_D grd_ph;
826
827 if (!fx) {
828 fx = space;
829 }
830
831 /* pressure stuff */
832 if (!qfast ||
833 qfast->bas_fcts != p_fe_space->bas_fcts || qfast->quad != quad) {
834 /* Aquire a QUAD_FAST structure for the evaluation of the gradient
835 * of the discrete pressure.
836 */
837 qfast = get_quad_fast(p_fe_space->bas_fcts, quad, INIT_GRD_PHI);
838
839 if (ph_loc) {
840 free_el_real_vec(ph_loc);
841 }
842 ph_loc = get_el_real_vec(p_fe_space->bas_fcts);
843 }
844
845 if (el != el_info->el) {
846 el = el_info->el;
847 fill_el_real_vec(ph_loc, el, p_h);
848 if (!(el_info->fill_flag & FILL_COORDS)) {
849 qelc = fill_quad_el_cache(el_info, quad,
850 FILL_EL_QUAD_WORLD|FILL_EL_QUAD_LAMBDA);
851 } else {
852 elgc = fill_el_geom_cache(el_info, FILL_EL_LAMBDA);
853 qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);
854 }
855 }
856
857 f(qelc->world[iq], fx);
858 SCAL_DOW(-1.0, fx);
859
860 /* 0-th order term */
861 AXPY_DOW(mu, uh_at_qp, fx);
862
863 if (!(el_info->fill_flag & FILL_COORDS)) {
864 eval_grd_uh_fast(grd_ph,
865 (const REAL_D *)qelc->param.Lambda[iq], ph_loc, qfast, iq);
866 } else {
867 eval_grd_uh_fast(grd_ph,
868 (const REAL_D *)elgc->Lambda, ph_loc, qfast, iq);
869 }
870
871 AXPY_DOW(1.0, grd_ph, fx);
872
873 return (const REAL *)fx;
874 }
875
876 /* Compute the residual for the stress boundary conditions.
877 *
878 * We need to compute (\nu D(u) - p I) n - g_{stress}.
879 */
stress_res(REAL_D result,const EL_INFO * el_info,const QUAD * quad,int qp,const REAL_D uh_qp,const REAL_D normal)880 static const REAL *stress_res(REAL_D result,
881 const EL_INFO *el_info,
882 const QUAD *quad, int qp,
883 const REAL_D uh_qp,
884 const REAL_D normal)
885 {
886 static REAL_D space;
887 static const QUAD_FAST *qfast;
888 static EL *el;
889 static EL_REAL_VEC *ph_loc;
890 REAL_D x;
891 REAL p_val;
892
893 if (!result) {
894 result = space;
895 }
896
897 /* pressure stuff */
898 if (!qfast || qfast->quad != quad ||
899 qfast->bas_fcts != p_fe_space->bas_fcts) {
900 /* Aquire a QUAD_FAST structure for the evaluation of the gradient
901 * of the discrete pressure.
902 */
903 qfast = get_quad_fast(p_fe_space->bas_fcts, quad, INIT_PHI);
904 if (ph_loc) {
905 free_el_real_vec(ph_loc);
906 }
907 ph_loc = get_el_real_vec(p_fe_space->bas_fcts);
908 }
909
910 if (el != el_info->el) {
911 el = el_info->el;
912 fill_el_real_vec(ph_loc, el, p_h);
913 INIT_ELEMENT(el_info, qfast);
914 }
915
916 p_val = eval_uh_fast(ph_loc, qfast, qp);
917
918 coord_to_world(el_info, quad->lambda[qp], x);
919
920 g_stress(x, normal, result);
921 AXPY_DOW(p_val, normal, result); /* "+" is correct here, the
922 * estimator only computes
923 * (A\nabla u_h + (A\nabla u_h)^t)
924 */
925
926 return result;
927 }
928
929 #define EOC(e, eo) log(eo/MAX(e, 1.0e-15))/M_LN2
930
estimate(MESH * mesh,ADAPT_STAT * adapt)931 static REAL estimate(MESH *mesh, ADAPT_STAT *adapt)
932 {
933 FUNCNAME("estimate");
934 static int init_done, only_dirichlet;
935 static REAL C[4] = { 1.0, 1.0, 0.0, 1.0 };
936 static REAL est_old = -1.0, u_err_old = -1.0, p_err_old = -1.0;
937 static REAL_DD A;
938 static const QUAD *quad;
939 static REAL est, u_err, p_err;
940 FLAGS fill_flags;
941 const void *est_handle;
942
943 if (!init_done) {
944 init_done = true;
945
946 /* Determine whether the pressure is determined only up to a
947 * constant.
948 */
949 only_dirichlet = true;
950 FOR_ALL_DOFS(bound->fe_space->admin,
951 if (bound->vec[dof] == NEUMANN) {
952 only_dirichlet = false;
953 }
954 );
955
956 GET_PARAMETER(1, "estimator C0", "%f", &C[0]);
957 GET_PARAMETER(1, "estimator C1", "%f", &C[1]);
958 GET_PARAMETER(1, "estimator C2", "%f", &C[2]);
959 GET_PARAMETER(1, "estimator C3", "%f", &C[3]);
960
961 quad = get_quadrature(mesh->dim, 2*u_fe_space->bas_fcts->degree);
962
963 A[0][0] = nu;
964 }
965
966 /* We use the standard error estimator for elliptic problems, but
967 * optionally add the L2-norm of the divergence of u_h, controlled
968 * by C[3].
969 */
970
971 /*********************** start of modified estimator ***********************/
972 est_handle = ellipt_est_dow_init(u_h, adapt, rw_el_est, NULL /* rw_estc */,
973 quad, NULL /* wall_quad */,
974 H1_NORM, C,
975 A, MATENT_REAL, MATENT_REAL,
976 true /* sym_grad */,
977 dirichlet_bndry,
978 r, INIT_GRD_UH|(mu != 0.0 ? INIT_UH : 0),
979 stress_res, 0);
980
981 fill_flags = FILL_NEIGH|FILL_COORDS|FILL_OPP_COORDS|FILL_BOUND|CALL_LEAF_EL;
982 TRAVERSE_FIRST(mesh, -1, fill_flags) {
983 const EL_GEOM_CACHE *elgc;
984 const QUAD_EL_CACHE *qelc;
985 REAL est_el;
986
987 est_el = element_est_d(el_info, est_handle);
988
989 if (C[3]) {
990 REAL div_uh_el, div_uh_qp;
991 const REAL_DD *grd_uh_qp;
992 int qp, i;
993
994 grd_uh_qp = element_est_grd_uh_d(est_handle);
995 div_uh_el = 0.0;
996 if (!(el_info->fill_flag & FILL_COORDS)) {
997 qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_DET);
998
999 for (qp = 0; qp < quad->n_points; qp++) {
1000 div_uh_qp = 0;
1001 for (i = 0; i < DIM_OF_WORLD; i++) {
1002 div_uh_qp += grd_uh_qp[qp][i][i];
1003 }
1004 div_uh_el += qelc->param.det[qp]*quad->w[qp]*SQR(div_uh_qp);
1005 }
1006 } else {
1007 elgc = fill_el_geom_cache(el_info, FILL_EL_DET);
1008
1009 for (qp = 0; qp < quad->n_points; qp++) {
1010 div_uh_qp = 0;
1011 for (i = 0; i < DIM_OF_WORLD; i++) {
1012 div_uh_qp += grd_uh_qp[qp][i][i];
1013 }
1014 div_uh_el += quad->w[qp]*SQR(div_uh_qp);
1015 }
1016 div_uh_el *= elgc->det;
1017 }
1018
1019 est_el += C[3] /* h2 */ * div_uh_el;
1020 }
1021
1022 element_est_d_finish(el_info, est_el, est_handle);
1023 } TRAVERSE_NEXT();
1024 est = ellipt_est_d_finish(adapt, est_handle);
1025 /*********************** end of modified estimator ***********************/
1026
1027 MSG("estimate = %.8le", est);
1028 if (est_old >= 0) {
1029 print_msg(", EOC: %.2lf\n", EOC(est,est_old));
1030 } else {
1031 print_msg("\n");
1032 }
1033 est_old = est;
1034
1035 quad = get_quadrature(mesh->dim, 2*(u_h->fe_space->bas_fcts->degree - 1));
1036
1037 u_err = H1_err_dow(grd_u, u_h, quad, 0, NULL, NULL);
1038
1039 MSG("||u-uh||H1 = %.8le", u_err);
1040 if (u_err_old >= 0) {
1041 print_msg(", EOC: %.2lf\n", EOC(u_err, u_err_old));
1042 } else {
1043 print_msg("\n");
1044 }
1045 u_err_old = u_err;
1046
1047 p_err = L2_err(p, p_h, quad, false /* rel_error */,
1048 only_dirichlet /* mean-value*/,
1049 NULL, NULL);
1050
1051 MSG("||p-ph||L2 = %.8le", p_err);
1052 if (p_err_old >= 0) {
1053 print_msg(", EOC: %.2lf\n", EOC(p_err, p_err_old));
1054 } else {
1055 print_msg("\n");
1056 }
1057 p_err_old = p_err;
1058
1059 MSG("(||u-uh||H1 + ||p-p_h||L2)/estimate = %.2lf\n",
1060 (u_err + p_err)/MAX(est,1.e-15));
1061
1062 if (do_graphics) {
1063 graphics_d(mesh, u_h, p_h, get_el_est, u, HUGE_VAL /* time */);
1064 }
1065
1066 return adapt->err_sum;
1067 }
1068
1069 /*******************************************************************************
1070 * The main program.
1071 ******************************************************************************/
1072
main(int argc,char ** argv)1073 int main(int argc, char **argv)
1074 {
1075 FUNCNAME("main");
1076 MACRO_DATA *data;
1077 MESH *mesh;
1078 int n_refine = 0, dim, u_degree = 2;
1079 int dirichlet_bit = 0;
1080 ADAPT_STAT *adapt;
1081 STOKES_PAIR stokes;
1082 char stokes_name[256];
1083 char filename[PATH_MAX];
1084
1085 /*****************************************************************************
1086 * first of all, initialize the access to parameters of the init file
1087 ****************************************************************************/
1088 parse_parameters(argc, argv, "INIT/quasi-stokes.dat");
1089
1090 /*****************************************************************************
1091 * then read some basic parameters
1092 ****************************************************************************/
1093 GET_PARAMETER(1, "mesh dimension", "%d", &dim);
1094 GET_PARAMETER(1, "macro file name", "%s", filename);
1095 GET_PARAMETER(1, "global refinements", "%d", &n_refine);
1096 GET_PARAMETER(1, "online graphics", "%d", &do_graphics);
1097 GET_PARAMETER(1, "dirichlet boundary", "%d", &dirichlet_bit);
1098
1099 GET_PARAMETER(1, "stokes pair", "%s", stokes_name);
1100 GET_PARAMETER(1, "velocity degree", "%d", &u_degree);
1101
1102 GET_PARAMETER(1, "QS mu", "%f", &mu);
1103 GET_PARAMETER(1, "QS nu", "%f", &nu);
1104
1105 GET_PARAMETER(1, "velocity translation",
1106 SCAN_FORMAT_DOW, SCAN_EXPAND_DOW(v_shift));
1107 GET_PARAMETER(1, "pressure translation",
1108 SCAN_FORMAT_DOW, SCAN_EXPAND_DOW(p_shift));
1109
1110 alpha_robin = mu; /* black magic, should probably be weighted by
1111 * some power of h
1112 */
1113
1114 /*****************************************************************************
1115 * get a mesh, and read the macro triangulation from file
1116 ****************************************************************************/
1117 data = read_macro(filename);
1118 mesh = GET_MESH(dim, "ALBERTA mesh", data, NULL, NULL);
1119 free_macro_data(data);
1120 init_leaf_data(mesh, sizeof(struct ellipt_leaf_data), NULL, NULL);
1121
1122 /*****************************************************************************
1123 * aquire velocity and pressure space finite element spaces
1124 ****************************************************************************/
1125 stokes = stokes_pair(stokes_name, dim, u_degree);
1126 u_fe_space = get_fe_space(mesh, "velocity", stokes.velocity, DIM_OF_WORLD,
1127 ADM_FLAGS_DFLT);
1128 p_fe_space = get_fe_space(mesh, "pressure", stokes.pressure, 1,
1129 ADM_FLAGS_DFLT);
1130
1131 /*****************************************************************************
1132 * Globally refine the mesh a little bit. Maybe do some graphics
1133 * output to amuse the spectators.
1134 ****************************************************************************/
1135 global_refine(mesh, n_refine * mesh->dim, FILL_NOTHING);
1136
1137 if (do_graphics) {
1138 graphics_d(mesh, NULL, NULL, NULL, NULL, HUGE_VAL /* time */);
1139 }
1140
1141 /*****************************************************************************
1142 * initialize the global variables shared across build(), solve()
1143 * and estimate().
1144 ****************************************************************************/
1145 A = get_dof_matrix("A", u_fe_space, NULL);
1146 dof_matrix_try_diagonal(A); /* partly diagonalize, e.g. for Mini */
1147 f_h = get_dof_real_vec_d("f_h", u_fe_space);
1148 u_h = get_dof_real_vec_d("u_h", u_fe_space);
1149 set_refine_inter_dow(u_h);
1150 set_coarse_inter_dow(u_h);
1151
1152 B = get_dof_matrix("B", u_fe_space, p_fe_space);
1153 dof_matrix_try_diagonal(B);
1154 Bt = get_dof_matrix("Bt", p_fe_space, u_fe_space);
1155 dof_matrix_try_diagonal(Bt);
1156 p_h = get_dof_real_vec("p_h", p_fe_space);
1157 set_refine_inter(p_h);
1158 set_coarse_inter(p_h);
1159
1160 g_h = get_dof_real_vec("g_h", p_fe_space);
1161 bound = get_dof_schar_vec("bound", u_fe_space);
1162
1163 Yproj = get_dof_matrix("Yproj", p_fe_space, NULL);
1164 dof_matrix_try_diagonal(Yproj);
1165 Yprec = get_dof_matrix("Yprec", p_fe_space, NULL);
1166 // Nope, Yprec is never structurally diagonal
1167
1168 /* Do not forget to initialize u_h and p_h */
1169 dof_set_dow(0.0, u_h);
1170 dof_set(0.0, p_h);
1171
1172 /* Boundary classifiction */
1173 BNDRY_FLAGS_ALL(stress_bndry);
1174 if (dirichlet_bit > 0) {
1175 BNDRY_FLAGS_SET(dirichlet_bndry, dirichlet_bit);
1176 BNDRY_FLAGS_XOR(stress_bndry, dirichlet_bndry);
1177 BNDRY_FLAGS_MARK_BNDRY(stress_bndry);
1178 }
1179
1180 /*****************************************************************************
1181 * init the adapt structure and start adaptive method
1182 ****************************************************************************/
1183 adapt = get_adapt_stat(mesh->dim, "Quasi-Stokes", "adapt", 2, NULL);
1184 adapt->estimate = estimate;
1185 adapt->get_el_est = get_el_est;
1186 adapt->build_after_coarsen = build;
1187 adapt->solve = solve;
1188
1189 adapt_method_stat(mesh, adapt);
1190
1191 WAIT_REALLY;
1192
1193 return 0;
1194 }
1195