Actual source code: ex42.c

      1: static char help[] = "Solves the incompressible, variable viscosity Stokes equation in 3d using Q1Q1 elements, \n\
      2: stabilized with Bochev's polynomial projection method. Note that implementation here assumes \n\
      3: all boundaries are free-slip, i.e. zero normal flow and zero tangential stress \n\
      4:      -mx : number elements in x-direction \n\
      5:      -my : number elements in y-direction \n\
      6:      -mz : number elements in z-direction \n\
      7:      -stokes_ksp_monitor_blocks : active monitor for each component u,v,w,p \n\
      8:      -model : defines viscosity and forcing function. Choose either: 0 (isoviscous), 1 (manufactured-broken), 2 (solcx-2d), 3 (sinker), 4 (subdomain jumps) \n";

     10: /* Contributed by Dave May */

     12: #include <petscksp.h>
     13: #include <petscdmda.h>

     15: #define PROFILE_TIMING
     16: #define ASSEMBLE_LOWER_TRIANGULAR

     18: #define NSD          3 /* number of spatial dimensions */
     19: #define NODES_PER_EL 8 /* nodes per element */
     20: #define U_DOFS       3 /* degrees of freedom per velocity node */
     21: #define P_DOFS       1 /* degrees of freedom per pressure node */
     22: #define GAUSS_POINTS 8

     24: /* Gauss point based evaluation */
     25: typedef struct {
     26:   PetscScalar gp_coords[NSD * GAUSS_POINTS];
     27:   PetscScalar eta[GAUSS_POINTS];
     28:   PetscScalar fx[GAUSS_POINTS];
     29:   PetscScalar fy[GAUSS_POINTS];
     30:   PetscScalar fz[GAUSS_POINTS];
     31:   PetscScalar hc[GAUSS_POINTS];
     32: } GaussPointCoefficients;

     34: typedef struct {
     35:   PetscScalar u_dof;
     36:   PetscScalar v_dof;
     37:   PetscScalar w_dof;
     38:   PetscScalar p_dof;
     39: } StokesDOF;

     41: typedef struct _p_CellProperties *CellProperties;
     42: struct _p_CellProperties {
     43:   PetscInt                ncells;
     44:   PetscInt                mx, my, mz;
     45:   PetscInt                sex, sey, sez;
     46:   GaussPointCoefficients *gpc;
     47: };

     49: /* elements */
     50: PetscErrorCode CellPropertiesCreate(DM da_stokes, CellProperties *C)
     51: {
     52:   CellProperties cells;
     53:   PetscInt       mx, my, mz, sex, sey, sez;

     55:   PetscFunctionBeginUser;
     56:   PetscCall(PetscNew(&cells));

     58:   PetscCall(DMDAGetElementsCorners(da_stokes, &sex, &sey, &sez));
     59:   PetscCall(DMDAGetElementsSizes(da_stokes, &mx, &my, &mz));
     60:   cells->mx     = mx;
     61:   cells->my     = my;
     62:   cells->mz     = mz;
     63:   cells->ncells = mx * my * mz;
     64:   cells->sex    = sex;
     65:   cells->sey    = sey;
     66:   cells->sez    = sez;

     68:   PetscCall(PetscMalloc1(mx * my * mz, &cells->gpc));

     70:   *C = cells;
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
     75: {
     76:   CellProperties cells;

     78:   PetscFunctionBeginUser;
     79:   if (!C) PetscFunctionReturn(PETSC_SUCCESS);
     80:   cells = *C;
     81:   PetscCall(PetscFree(cells->gpc));
     82:   PetscCall(PetscFree(cells));
     83:   *C = NULL;
     84:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

     87: PetscErrorCode CellPropertiesGetCell(CellProperties C, PetscInt II, PetscInt J, PetscInt K, GaussPointCoefficients **G)
     88: {
     89:   PetscFunctionBeginUser;
     90:   *G = &C->gpc[(II - C->sex) + (J - C->sey) * C->mx + (K - C->sez) * C->mx * C->my];
     91:   PetscFunctionReturn(PETSC_SUCCESS);
     92: }

     94: /* FEM routines */
     95: /*
     96:  Element: Local basis function ordering
     97:  1-----2
     98:  |     |
     99:  |     |
    100:  0-----3
    101:  */
    102: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[], PetscScalar Ni[])
    103: {
    104:   PetscReal xi   = PetscRealPart(_xi[0]);
    105:   PetscReal eta  = PetscRealPart(_xi[1]);
    106:   PetscReal zeta = PetscRealPart(_xi[2]);

    108:   Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
    109:   Ni[1] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
    110:   Ni[2] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
    111:   Ni[3] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);

    113:   Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
    114:   Ni[5] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
    115:   Ni[6] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
    116:   Ni[7] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
    117: }

    119: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
    120: {
    121:   PetscReal xi   = PetscRealPart(_xi[0]);
    122:   PetscReal eta  = PetscRealPart(_xi[1]);
    123:   PetscReal zeta = PetscRealPart(_xi[2]);
    124:   /* xi deriv */
    125:   GNi[0][0] = -0.125 * (1.0 - eta) * (1.0 - zeta);
    126:   GNi[0][1] = -0.125 * (1.0 + eta) * (1.0 - zeta);
    127:   GNi[0][2] = 0.125 * (1.0 + eta) * (1.0 - zeta);
    128:   GNi[0][3] = 0.125 * (1.0 - eta) * (1.0 - zeta);

    130:   GNi[0][4] = -0.125 * (1.0 - eta) * (1.0 + zeta);
    131:   GNi[0][5] = -0.125 * (1.0 + eta) * (1.0 + zeta);
    132:   GNi[0][6] = 0.125 * (1.0 + eta) * (1.0 + zeta);
    133:   GNi[0][7] = 0.125 * (1.0 - eta) * (1.0 + zeta);
    134:   /* eta deriv */
    135:   GNi[1][0] = -0.125 * (1.0 - xi) * (1.0 - zeta);
    136:   GNi[1][1] = 0.125 * (1.0 - xi) * (1.0 - zeta);
    137:   GNi[1][2] = 0.125 * (1.0 + xi) * (1.0 - zeta);
    138:   GNi[1][3] = -0.125 * (1.0 + xi) * (1.0 - zeta);

    140:   GNi[1][4] = -0.125 * (1.0 - xi) * (1.0 + zeta);
    141:   GNi[1][5] = 0.125 * (1.0 - xi) * (1.0 + zeta);
    142:   GNi[1][6] = 0.125 * (1.0 + xi) * (1.0 + zeta);
    143:   GNi[1][7] = -0.125 * (1.0 + xi) * (1.0 + zeta);
    144:   /* zeta deriv */
    145:   GNi[2][0] = -0.125 * (1.0 - xi) * (1.0 - eta);
    146:   GNi[2][1] = -0.125 * (1.0 - xi) * (1.0 + eta);
    147:   GNi[2][2] = -0.125 * (1.0 + xi) * (1.0 + eta);
    148:   GNi[2][3] = -0.125 * (1.0 + xi) * (1.0 - eta);

    150:   GNi[2][4] = 0.125 * (1.0 - xi) * (1.0 - eta);
    151:   GNi[2][5] = 0.125 * (1.0 - xi) * (1.0 + eta);
    152:   GNi[2][6] = 0.125 * (1.0 + xi) * (1.0 + eta);
    153:   GNi[2][7] = 0.125 * (1.0 + xi) * (1.0 - eta);
    154: }

    156: static void matrix_inverse_3x3(PetscScalar A[3][3], PetscScalar B[3][3])
    157: {
    158:   PetscScalar t4, t6, t8, t10, t12, t14, t17;

    160:   t4  = A[2][0] * A[0][1];
    161:   t6  = A[2][0] * A[0][2];
    162:   t8  = A[1][0] * A[0][1];
    163:   t10 = A[1][0] * A[0][2];
    164:   t12 = A[0][0] * A[1][1];
    165:   t14 = A[0][0] * A[1][2];
    166:   t17 = 0.1e1 / (t4 * A[1][2] - t6 * A[1][1] - t8 * A[2][2] + t10 * A[2][1] + t12 * A[2][2] - t14 * A[2][1]);

    168:   B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
    169:   B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
    170:   B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
    171:   B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
    172:   B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
    173:   B[1][2] = -(-t10 + t14) * t17;
    174:   B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
    175:   B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
    176:   B[2][2] = (-t8 + t12) * t17;
    177: }

    179: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
    180: {
    181:   PetscScalar J00, J01, J02, J10, J11, J12, J20, J21, J22;
    182:   PetscInt    n;
    183:   PetscScalar iJ[3][3], JJ[3][3];

    185:   J00 = J01 = J02 = 0.0;
    186:   J10 = J11 = J12 = 0.0;
    187:   J20 = J21 = J22 = 0.0;
    188:   for (n = 0; n < NODES_PER_EL; n++) {
    189:     PetscScalar cx = coords[NSD * n + 0];
    190:     PetscScalar cy = coords[NSD * n + 1];
    191:     PetscScalar cz = coords[NSD * n + 2];

    193:     /* J_ij = d(x_j) / d(xi_i) */ /* J_ij = \sum _I GNi[j][I} * x_i */
    194:     J00 = J00 + GNi[0][n] * cx;   /* J_xx */
    195:     J01 = J01 + GNi[0][n] * cy;   /* J_xy = dx/deta */
    196:     J02 = J02 + GNi[0][n] * cz;   /* J_xz = dx/dzeta */

    198:     J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
    199:     J11 = J11 + GNi[1][n] * cy; /* J_yy */
    200:     J12 = J12 + GNi[1][n] * cz; /* J_yz */

    202:     J20 = J20 + GNi[2][n] * cx; /* J_zx */
    203:     J21 = J21 + GNi[2][n] * cy; /* J_zy */
    204:     J22 = J22 + GNi[2][n] * cz; /* J_zz */
    205:   }

    207:   JJ[0][0] = J00;
    208:   JJ[0][1] = J01;
    209:   JJ[0][2] = J02;
    210:   JJ[1][0] = J10;
    211:   JJ[1][1] = J11;
    212:   JJ[1][2] = J12;
    213:   JJ[2][0] = J20;
    214:   JJ[2][1] = J21;
    215:   JJ[2][2] = J22;

    217:   matrix_inverse_3x3(JJ, iJ);

    219:   *det_J = J00 * J11 * J22 - J00 * J12 * J21 - J10 * J01 * J22 + J10 * J02 * J21 + J20 * J01 * J12 - J20 * J02 * J11;

    221:   for (n = 0; n < NODES_PER_EL; n++) {
    222:     GNx[0][n] = GNi[0][n] * iJ[0][0] + GNi[1][n] * iJ[0][1] + GNi[2][n] * iJ[0][2];
    223:     GNx[1][n] = GNi[0][n] * iJ[1][0] + GNi[1][n] * iJ[1][1] + GNi[2][n] * iJ[1][2];
    224:     GNx[2][n] = GNi[0][n] * iJ[2][0] + GNi[1][n] * iJ[2][1] + GNi[2][n] * iJ[2][2];
    225:   }
    226: }

    228: static void ConstructGaussQuadrature3D(PetscInt *ngp, PetscScalar gp_xi[][NSD], PetscScalar gp_weight[])
    229: {
    230:   *ngp        = 8;
    231:   gp_xi[0][0] = -0.57735026919;
    232:   gp_xi[0][1] = -0.57735026919;
    233:   gp_xi[0][2] = -0.57735026919;
    234:   gp_xi[1][0] = -0.57735026919;
    235:   gp_xi[1][1] = 0.57735026919;
    236:   gp_xi[1][2] = -0.57735026919;
    237:   gp_xi[2][0] = 0.57735026919;
    238:   gp_xi[2][1] = 0.57735026919;
    239:   gp_xi[2][2] = -0.57735026919;
    240:   gp_xi[3][0] = 0.57735026919;
    241:   gp_xi[3][1] = -0.57735026919;
    242:   gp_xi[3][2] = -0.57735026919;

    244:   gp_xi[4][0] = -0.57735026919;
    245:   gp_xi[4][1] = -0.57735026919;
    246:   gp_xi[4][2] = 0.57735026919;
    247:   gp_xi[5][0] = -0.57735026919;
    248:   gp_xi[5][1] = 0.57735026919;
    249:   gp_xi[5][2] = 0.57735026919;
    250:   gp_xi[6][0] = 0.57735026919;
    251:   gp_xi[6][1] = 0.57735026919;
    252:   gp_xi[6][2] = 0.57735026919;
    253:   gp_xi[7][0] = 0.57735026919;
    254:   gp_xi[7][1] = -0.57735026919;
    255:   gp_xi[7][2] = 0.57735026919;

    257:   gp_weight[0] = 1.0;
    258:   gp_weight[1] = 1.0;
    259:   gp_weight[2] = 1.0;
    260:   gp_weight[3] = 1.0;

    262:   gp_weight[4] = 1.0;
    263:   gp_weight[5] = 1.0;
    264:   gp_weight[6] = 1.0;
    265:   gp_weight[7] = 1.0;
    266: }

    268: /*
    269:  i,j are the element indices
    270:  The unknown is a vector quantity.
    271:  The s[].c is used to indicate the degree of freedom.
    272:  */
    273: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[], MatStencil s_p[], PetscInt i, PetscInt j, PetscInt k)
    274: {
    275:   PetscInt n;

    277:   PetscFunctionBeginUser;
    278:   /* velocity */
    279:   n = 0;
    280:   /* node 0 */
    281:   s_u[n].i = i;
    282:   s_u[n].j = j;
    283:   s_u[n].k = k;
    284:   s_u[n].c = 0;
    285:   n++; /* Vx0 */
    286:   s_u[n].i = i;
    287:   s_u[n].j = j;
    288:   s_u[n].k = k;
    289:   s_u[n].c = 1;
    290:   n++; /* Vy0 */
    291:   s_u[n].i = i;
    292:   s_u[n].j = j;
    293:   s_u[n].k = k;
    294:   s_u[n].c = 2;
    295:   n++; /* Vz0 */

    297:   s_u[n].i = i;
    298:   s_u[n].j = j + 1;
    299:   s_u[n].k = k;
    300:   s_u[n].c = 0;
    301:   n++;
    302:   s_u[n].i = i;
    303:   s_u[n].j = j + 1;
    304:   s_u[n].k = k;
    305:   s_u[n].c = 1;
    306:   n++;
    307:   s_u[n].i = i;
    308:   s_u[n].j = j + 1;
    309:   s_u[n].k = k;
    310:   s_u[n].c = 2;
    311:   n++;

    313:   s_u[n].i = i + 1;
    314:   s_u[n].j = j + 1;
    315:   s_u[n].k = k;
    316:   s_u[n].c = 0;
    317:   n++;
    318:   s_u[n].i = i + 1;
    319:   s_u[n].j = j + 1;
    320:   s_u[n].k = k;
    321:   s_u[n].c = 1;
    322:   n++;
    323:   s_u[n].i = i + 1;
    324:   s_u[n].j = j + 1;
    325:   s_u[n].k = k;
    326:   s_u[n].c = 2;
    327:   n++;

    329:   s_u[n].i = i + 1;
    330:   s_u[n].j = j;
    331:   s_u[n].k = k;
    332:   s_u[n].c = 0;
    333:   n++;
    334:   s_u[n].i = i + 1;
    335:   s_u[n].j = j;
    336:   s_u[n].k = k;
    337:   s_u[n].c = 1;
    338:   n++;
    339:   s_u[n].i = i + 1;
    340:   s_u[n].j = j;
    341:   s_u[n].k = k;
    342:   s_u[n].c = 2;
    343:   n++;

    345:   /* */
    346:   s_u[n].i = i;
    347:   s_u[n].j = j;
    348:   s_u[n].k = k + 1;
    349:   s_u[n].c = 0;
    350:   n++; /* Vx4 */
    351:   s_u[n].i = i;
    352:   s_u[n].j = j;
    353:   s_u[n].k = k + 1;
    354:   s_u[n].c = 1;
    355:   n++; /* Vy4 */
    356:   s_u[n].i = i;
    357:   s_u[n].j = j;
    358:   s_u[n].k = k + 1;
    359:   s_u[n].c = 2;
    360:   n++; /* Vz4 */

    362:   s_u[n].i = i;
    363:   s_u[n].j = j + 1;
    364:   s_u[n].k = k + 1;
    365:   s_u[n].c = 0;
    366:   n++;
    367:   s_u[n].i = i;
    368:   s_u[n].j = j + 1;
    369:   s_u[n].k = k + 1;
    370:   s_u[n].c = 1;
    371:   n++;
    372:   s_u[n].i = i;
    373:   s_u[n].j = j + 1;
    374:   s_u[n].k = k + 1;
    375:   s_u[n].c = 2;
    376:   n++;

    378:   s_u[n].i = i + 1;
    379:   s_u[n].j = j + 1;
    380:   s_u[n].k = k + 1;
    381:   s_u[n].c = 0;
    382:   n++;
    383:   s_u[n].i = i + 1;
    384:   s_u[n].j = j + 1;
    385:   s_u[n].k = k + 1;
    386:   s_u[n].c = 1;
    387:   n++;
    388:   s_u[n].i = i + 1;
    389:   s_u[n].j = j + 1;
    390:   s_u[n].k = k + 1;
    391:   s_u[n].c = 2;
    392:   n++;

    394:   s_u[n].i = i + 1;
    395:   s_u[n].j = j;
    396:   s_u[n].k = k + 1;
    397:   s_u[n].c = 0;
    398:   n++;
    399:   s_u[n].i = i + 1;
    400:   s_u[n].j = j;
    401:   s_u[n].k = k + 1;
    402:   s_u[n].c = 1;
    403:   n++;
    404:   s_u[n].i = i + 1;
    405:   s_u[n].j = j;
    406:   s_u[n].k = k + 1;
    407:   s_u[n].c = 2;
    408:   n++;

    410:   /* pressure */
    411:   n = 0;

    413:   s_p[n].i = i;
    414:   s_p[n].j = j;
    415:   s_p[n].k = k;
    416:   s_p[n].c = 3;
    417:   n++; /* P0 */
    418:   s_p[n].i = i;
    419:   s_p[n].j = j + 1;
    420:   s_p[n].k = k;
    421:   s_p[n].c = 3;
    422:   n++;
    423:   s_p[n].i = i + 1;
    424:   s_p[n].j = j + 1;
    425:   s_p[n].k = k;
    426:   s_p[n].c = 3;
    427:   n++;
    428:   s_p[n].i = i + 1;
    429:   s_p[n].j = j;
    430:   s_p[n].k = k;
    431:   s_p[n].c = 3;
    432:   n++;

    434:   s_p[n].i = i;
    435:   s_p[n].j = j;
    436:   s_p[n].k = k + 1;
    437:   s_p[n].c = 3;
    438:   n++; /* P0 */
    439:   s_p[n].i = i;
    440:   s_p[n].j = j + 1;
    441:   s_p[n].k = k + 1;
    442:   s_p[n].c = 3;
    443:   n++;
    444:   s_p[n].i = i + 1;
    445:   s_p[n].j = j + 1;
    446:   s_p[n].k = k + 1;
    447:   s_p[n].c = 3;
    448:   n++;
    449:   s_p[n].i = i + 1;
    450:   s_p[n].j = j;
    451:   s_p[n].k = k + 1;
    452:   s_p[n].c = 3;
    453:   n++;
    454:   PetscFunctionReturn(PETSC_SUCCESS);
    455: }

    457: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords, PetscInt i, PetscInt j, PetscInt k, PetscScalar el_coord[])
    458: {
    459:   PetscFunctionBeginUser;
    460:   /* get coords for the element */
    461:   el_coord[0] = coords[k][j][i].x;
    462:   el_coord[1] = coords[k][j][i].y;
    463:   el_coord[2] = coords[k][j][i].z;

    465:   el_coord[3] = coords[k][j + 1][i].x;
    466:   el_coord[4] = coords[k][j + 1][i].y;
    467:   el_coord[5] = coords[k][j + 1][i].z;

    469:   el_coord[6] = coords[k][j + 1][i + 1].x;
    470:   el_coord[7] = coords[k][j + 1][i + 1].y;
    471:   el_coord[8] = coords[k][j + 1][i + 1].z;

    473:   el_coord[9]  = coords[k][j][i + 1].x;
    474:   el_coord[10] = coords[k][j][i + 1].y;
    475:   el_coord[11] = coords[k][j][i + 1].z;

    477:   el_coord[12] = coords[k + 1][j][i].x;
    478:   el_coord[13] = coords[k + 1][j][i].y;
    479:   el_coord[14] = coords[k + 1][j][i].z;

    481:   el_coord[15] = coords[k + 1][j + 1][i].x;
    482:   el_coord[16] = coords[k + 1][j + 1][i].y;
    483:   el_coord[17] = coords[k + 1][j + 1][i].z;

    485:   el_coord[18] = coords[k + 1][j + 1][i + 1].x;
    486:   el_coord[19] = coords[k + 1][j + 1][i + 1].y;
    487:   el_coord[20] = coords[k + 1][j + 1][i + 1].z;

    489:   el_coord[21] = coords[k + 1][j][i + 1].x;
    490:   el_coord[22] = coords[k + 1][j][i + 1].y;
    491:   el_coord[23] = coords[k + 1][j][i + 1].z;
    492:   PetscFunctionReturn(PETSC_SUCCESS);
    493: }

    495: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field, PetscInt i, PetscInt j, PetscInt k, StokesDOF nodal_fields[])
    496: {
    497:   PetscFunctionBeginUser;
    498:   /* get the nodal fields for u */
    499:   nodal_fields[0].u_dof = field[k][j][i].u_dof;
    500:   nodal_fields[0].v_dof = field[k][j][i].v_dof;
    501:   nodal_fields[0].w_dof = field[k][j][i].w_dof;

    503:   nodal_fields[1].u_dof = field[k][j + 1][i].u_dof;
    504:   nodal_fields[1].v_dof = field[k][j + 1][i].v_dof;
    505:   nodal_fields[1].w_dof = field[k][j + 1][i].w_dof;

    507:   nodal_fields[2].u_dof = field[k][j + 1][i + 1].u_dof;
    508:   nodal_fields[2].v_dof = field[k][j + 1][i + 1].v_dof;
    509:   nodal_fields[2].w_dof = field[k][j + 1][i + 1].w_dof;

    511:   nodal_fields[3].u_dof = field[k][j][i + 1].u_dof;
    512:   nodal_fields[3].v_dof = field[k][j][i + 1].v_dof;
    513:   nodal_fields[3].w_dof = field[k][j][i + 1].w_dof;

    515:   nodal_fields[4].u_dof = field[k + 1][j][i].u_dof;
    516:   nodal_fields[4].v_dof = field[k + 1][j][i].v_dof;
    517:   nodal_fields[4].w_dof = field[k + 1][j][i].w_dof;

    519:   nodal_fields[5].u_dof = field[k + 1][j + 1][i].u_dof;
    520:   nodal_fields[5].v_dof = field[k + 1][j + 1][i].v_dof;
    521:   nodal_fields[5].w_dof = field[k + 1][j + 1][i].w_dof;

    523:   nodal_fields[6].u_dof = field[k + 1][j + 1][i + 1].u_dof;
    524:   nodal_fields[6].v_dof = field[k + 1][j + 1][i + 1].v_dof;
    525:   nodal_fields[6].w_dof = field[k + 1][j + 1][i + 1].w_dof;

    527:   nodal_fields[7].u_dof = field[k + 1][j][i + 1].u_dof;
    528:   nodal_fields[7].v_dof = field[k + 1][j][i + 1].v_dof;
    529:   nodal_fields[7].w_dof = field[k + 1][j][i + 1].w_dof;

    531:   /* get the nodal fields for p */
    532:   nodal_fields[0].p_dof = field[k][j][i].p_dof;
    533:   nodal_fields[1].p_dof = field[k][j + 1][i].p_dof;
    534:   nodal_fields[2].p_dof = field[k][j + 1][i + 1].p_dof;
    535:   nodal_fields[3].p_dof = field[k][j][i + 1].p_dof;

    537:   nodal_fields[4].p_dof = field[k + 1][j][i].p_dof;
    538:   nodal_fields[5].p_dof = field[k + 1][j + 1][i].p_dof;
    539:   nodal_fields[6].p_dof = field[k + 1][j + 1][i + 1].p_dof;
    540:   nodal_fields[7].p_dof = field[k + 1][j][i + 1].p_dof;
    541:   PetscFunctionReturn(PETSC_SUCCESS);
    542: }

    544: static PetscInt ASS_MAP_wIwDI_uJuDJ(PetscInt wi, PetscInt wd, PetscInt w_NPE, PetscInt w_dof, PetscInt ui, PetscInt ud, PetscInt u_NPE, PetscInt u_dof)
    545: {
    546:   PetscInt              ij;
    547:   PETSC_UNUSED PetscInt r, c, nr, nc;

    549:   nr = w_NPE * w_dof;
    550:   nc = u_NPE * u_dof;

    552:   r = w_dof * wi + wd;
    553:   c = u_dof * ui + ud;

    555:   ij = r * nc + c;

    557:   return ij;
    558: }

    560: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
    561: {
    562:   PetscInt n, II, J, K;

    564:   PetscFunctionBeginUser;
    565:   for (n = 0; n < NODES_PER_EL; n++) {
    566:     II = u_eqn[NSD * n].i;
    567:     J  = u_eqn[NSD * n].j;
    568:     K  = u_eqn[NSD * n].k;

    570:     fields_F[K][J][II].u_dof = fields_F[K][J][II].u_dof + Fe_u[NSD * n];

    572:     II = u_eqn[NSD * n + 1].i;
    573:     J  = u_eqn[NSD * n + 1].j;
    574:     K  = u_eqn[NSD * n + 1].k;

    576:     fields_F[K][J][II].v_dof = fields_F[K][J][II].v_dof + Fe_u[NSD * n + 1];

    578:     II                       = u_eqn[NSD * n + 2].i;
    579:     J                        = u_eqn[NSD * n + 2].j;
    580:     K                        = u_eqn[NSD * n + 2].k;
    581:     fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof + Fe_u[NSD * n + 2];

    583:     II = p_eqn[n].i;
    584:     J  = p_eqn[n].j;
    585:     K  = p_eqn[n].k;

    587:     fields_F[K][J][II].p_dof = fields_F[K][J][II].p_dof + Fe_p[n];
    588:   }
    589:   PetscFunctionReturn(PETSC_SUCCESS);
    590: }

    592: static void FormStressOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    593: {
    594:   PetscInt       ngp;
    595:   PetscScalar    gp_xi[GAUSS_POINTS][NSD];
    596:   PetscScalar    gp_weight[GAUSS_POINTS];
    597:   PetscInt       p, i, j, k;
    598:   PetscScalar    GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    599:   PetscScalar    J_p, tildeD[6];
    600:   PetscScalar    B[6][U_DOFS * NODES_PER_EL];
    601:   const PetscInt nvdof = U_DOFS * NODES_PER_EL;

    603:   /* define quadrature rule */
    604:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    606:   /* evaluate integral */
    607:   for (p = 0; p < ngp; p++) {
    608:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    609:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);

    611:     for (i = 0; i < NODES_PER_EL; i++) {
    612:       PetscScalar d_dx_i = GNx_p[0][i];
    613:       PetscScalar d_dy_i = GNx_p[1][i];
    614:       PetscScalar d_dz_i = GNx_p[2][i];

    616:       B[0][3 * i]     = d_dx_i;
    617:       B[0][3 * i + 1] = 0.0;
    618:       B[0][3 * i + 2] = 0.0;
    619:       B[1][3 * i]     = 0.0;
    620:       B[1][3 * i + 1] = d_dy_i;
    621:       B[1][3 * i + 2] = 0.0;
    622:       B[2][3 * i]     = 0.0;
    623:       B[2][3 * i + 1] = 0.0;
    624:       B[2][3 * i + 2] = d_dz_i;

    626:       B[3][3 * i]     = d_dy_i;
    627:       B[3][3 * i + 1] = d_dx_i;
    628:       B[3][3 * i + 2] = 0.0; /* e_xy */
    629:       B[4][3 * i]     = d_dz_i;
    630:       B[4][3 * i + 1] = 0.0;
    631:       B[4][3 * i + 2] = d_dx_i; /* e_xz */
    632:       B[5][3 * i]     = 0.0;
    633:       B[5][3 * i + 1] = d_dz_i;
    634:       B[5][3 * i + 2] = d_dy_i; /* e_yz */
    635:     }

    637:     tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
    638:     tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
    639:     tildeD[2] = 2.0 * gp_weight[p] * J_p * eta[p];

    641:     tildeD[3] = gp_weight[p] * J_p * eta[p];
    642:     tildeD[4] = gp_weight[p] * J_p * eta[p];
    643:     tildeD[5] = gp_weight[p] * J_p * eta[p];

    645:     /* form Bt tildeD B */
    646:     /*
    647:      Ke_ij = Bt_ik . D_kl . B_lj
    648:      = B_ki . D_kl . B_lj
    649:      = B_ki . D_kk . B_kj
    650:      */
    651:     for (i = 0; i < nvdof; i++) {
    652:       for (j = i; j < nvdof; j++) {
    653:         for (k = 0; k < 6; k++) Ke[i * nvdof + j] += B[k][i] * tildeD[k] * B[k][j];
    654:       }
    655:     }
    656:   }
    657:   /* fill lower triangular part */
    658: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
    659:   for (i = 0; i < nvdof; i++) {
    660:     for (j = i; j < nvdof; j++) Ke[j * nvdof + i] = Ke[i * nvdof + j];
    661:   }
    662: #endif
    663: }

    665: static void FormGradientOperatorQ13D(PetscScalar Ke[], PetscScalar coords[])
    666: {
    667:   PetscInt    ngp;
    668:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
    669:   PetscScalar gp_weight[GAUSS_POINTS];
    670:   PetscInt    p, i, j, di;
    671:   PetscScalar Ni_p[NODES_PER_EL];
    672:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    673:   PetscScalar J_p, fac;

    675:   /* define quadrature rule */
    676:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    678:   /* evaluate integral */
    679:   for (p = 0; p < ngp; p++) {
    680:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
    681:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    682:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
    683:     fac = gp_weight[p] * J_p;

    685:     for (i = 0; i < NODES_PER_EL; i++) {     /* u nodes */
    686:       for (di = 0; di < NSD; di++) {         /* u dofs */
    687:         for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
    688:           PetscInt IJ;
    689:           IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 3, j, 0, NODES_PER_EL, 1);

    691:           Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
    692:         }
    693:       }
    694:     }
    695:   }
    696: }

    698: static void FormDivergenceOperatorQ13D(PetscScalar De[], PetscScalar coords[])
    699: {
    700:   PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
    701:   PetscInt    i, j;
    702:   PetscInt    nr_g, nc_g;

    704:   PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
    705:   FormGradientOperatorQ13D(Ge, coords);

    707:   nr_g = U_DOFS * NODES_PER_EL;
    708:   nc_g = P_DOFS * NODES_PER_EL;

    710:   for (i = 0; i < nr_g; i++) {
    711:     for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
    712:   }
    713: }

    715: static void FormStabilisationOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    716: {
    717:   PetscInt    ngp;
    718:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
    719:   PetscScalar gp_weight[GAUSS_POINTS];
    720:   PetscInt    p, i, j;
    721:   PetscScalar Ni_p[NODES_PER_EL];
    722:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    723:   PetscScalar J_p, fac, eta_avg;

    725:   /* define quadrature rule */
    726:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    728:   /* evaluate integral */
    729:   for (p = 0; p < ngp; p++) {
    730:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
    731:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    732:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
    733:     fac = gp_weight[p] * J_p;
    734:     /*
    735:      for (i = 0; i < NODES_PER_EL; i++) {
    736:      for (j = i; j < NODES_PER_EL; j++) {
    737:      Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
    738:      }
    739:      }
    740:      */

    742:     for (i = 0; i < NODES_PER_EL; i++) {
    743:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.015625);
    744:     }
    745:   }

    747:   /* scale */
    748:   eta_avg = 0.0;
    749:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    750:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
    751:   fac     = 1.0 / eta_avg;

    753:   /*
    754:    for (i = 0; i < NODES_PER_EL; i++) {
    755:    for (j = i; j < NODES_PER_EL; j++) {
    756:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
    757:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
    758:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
    759:    #endif
    760:    }
    761:    }
    762:    */

    764:   for (i = 0; i < NODES_PER_EL; i++) {
    765:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
    766:   }
    767: }

    769: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    770: {
    771:   PetscInt    ngp;
    772:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
    773:   PetscScalar gp_weight[GAUSS_POINTS];
    774:   PetscInt    p, i, j;
    775:   PetscScalar Ni_p[NODES_PER_EL];
    776:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    777:   PetscScalar J_p, fac, eta_avg;

    779:   /* define quadrature rule */
    780:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    782:   /* evaluate integral */
    783:   for (p = 0; p < ngp; p++) {
    784:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
    785:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    786:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
    787:     fac = gp_weight[p] * J_p;

    789:     /*
    790:      for (i = 0; i < NODES_PER_EL; i++) {
    791:      for (j = i; j < NODES_PER_EL; j++) {
    792:      Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
    793:      }
    794:      }
    795:      */

    797:     for (i = 0; i < NODES_PER_EL; i++) {
    798:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * Ni_p[i] * Ni_p[j];
    799:     }
    800:   }

    802:   /* scale */
    803:   eta_avg = 0.0;
    804:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    805:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
    806:   fac     = 1.0 / eta_avg;
    807:   /*
    808:    for (i = 0; i < NODES_PER_EL; i++) {
    809:    for (j = i; j < NODES_PER_EL; j++) {
    810:    Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
    811:    #if defined(ASSEMBLE_LOWER_TRIANGULAR)
    812:    Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
    813:    #endif
    814:    }
    815:    }
    816:    */

    818:   for (i = 0; i < NODES_PER_EL; i++) {
    819:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
    820:   }
    821: }

    823: static void FormMomentumRhsQ13D(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[], PetscScalar fz[])
    824: {
    825:   PetscInt    ngp;
    826:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
    827:   PetscScalar gp_weight[GAUSS_POINTS];
    828:   PetscInt    p, i;
    829:   PetscScalar Ni_p[NODES_PER_EL];
    830:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    831:   PetscScalar J_p, fac;

    833:   /* define quadrature rule */
    834:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    836:   /* evaluate integral */
    837:   for (p = 0; p < ngp; p++) {
    838:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
    839:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    840:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
    841:     fac = gp_weight[p] * J_p;

    843:     for (i = 0; i < NODES_PER_EL; i++) {
    844:       Fe[NSD * i] -= fac * Ni_p[i] * fx[p];
    845:       Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
    846:       Fe[NSD * i + 2] -= fac * Ni_p[i] * fz[p];
    847:     }
    848:   }
    849: }

    851: static void FormContinuityRhsQ13D(PetscScalar Fe[], PetscScalar coords[], PetscScalar hc[])
    852: {
    853:   PetscInt    ngp;
    854:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
    855:   PetscScalar gp_weight[GAUSS_POINTS];
    856:   PetscInt    p, i;
    857:   PetscScalar Ni_p[NODES_PER_EL];
    858:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    859:   PetscScalar J_p, fac;

    861:   /* define quadrature rule */
    862:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    864:   /* evaluate integral */
    865:   for (p = 0; p < ngp; p++) {
    866:     ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
    867:     ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    868:     ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
    869:     fac = gp_weight[p] * J_p;

    871:     for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac * Ni_p[i] * hc[p];
    872:   }
    873: }

    875: #define _ZERO_ROWCOL_i(A, i) \
    876:   do { \
    877:     PetscInt    KK; \
    878:     PetscScalar tmp = A[24 * (i) + (i)]; \
    879:     for (KK = 0; KK < 24; KK++) A[24 * (i) + KK] = 0.0; \
    880:     for (KK = 0; KK < 24; KK++) A[24 * KK + (i)] = 0.0; \
    881:     A[24 * (i) + (i)] = tmp; \
    882:   } while (0)

    884: #define _ZERO_ROW_i(A, i) \
    885:   do { \
    886:     PetscInt KK; \
    887:     for (KK = 0; KK < 8; KK++) A[8 * (i) + KK] = 0.0; \
    888:   } while (0)

    890: #define _ZERO_COL_i(A, i) \
    891:   do { \
    892:     PetscInt KK; \
    893:     for (KK = 0; KK < 8; KK++) A[24 * KK + (i)] = 0.0; \
    894:   } while (0)

    896: static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, CellProperties cell_properties)
    897: {
    898:   DM                      cda;
    899:   Vec                     coords;
    900:   DMDACoor3d           ***_coords;
    901:   MatStencil              u_eqn[NODES_PER_EL * U_DOFS];
    902:   MatStencil              p_eqn[NODES_PER_EL * P_DOFS];
    903:   PetscInt                sex, sey, sez, mx, my, mz;
    904:   PetscInt                ei, ej, ek;
    905:   PetscScalar             Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
    906:   PetscScalar             Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
    907:   PetscScalar             De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
    908:   PetscScalar             Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
    909:   PetscScalar             el_coords[NODES_PER_EL * NSD];
    910:   GaussPointCoefficients *props;
    911:   PetscScalar            *prop_eta;
    912:   PetscInt                n, M, N, P;

    914:   PetscFunctionBeginUser;
    915:   PetscCall(DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    916:   /* setup for coords */
    917:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    918:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    919:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    920:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
    921:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
    922:   for (ek = sez; ek < sez + mz; ek++) {
    923:     for (ej = sey; ej < sey + my; ej++) {
    924:       for (ei = sex; ei < sex + mx; ei++) {
    925:         /* get coords for the element */
    926:         PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
    927:         /* get cell properties */
    928:         PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &props));
    929:         /* get coefficients for the element */
    930:         prop_eta = props->eta;

    932:         /* initialise element stiffness matrix */
    933:         PetscCall(PetscMemzero(Ae, sizeof(Ae)));
    934:         PetscCall(PetscMemzero(Ge, sizeof(Ge)));
    935:         PetscCall(PetscMemzero(De, sizeof(De)));
    936:         PetscCall(PetscMemzero(Ce, sizeof(Ce)));

    938:         /* form element stiffness matrix */
    939:         FormStressOperatorQ13D(Ae, el_coords, prop_eta);
    940:         FormGradientOperatorQ13D(Ge, el_coords);
    941:         /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
    942:         FormDivergenceOperatorQ13D(De, el_coords);
    943:         /*#endif*/
    944:         FormStabilisationOperatorQ13D(Ce, el_coords, prop_eta);

    946:         /* insert element matrix into global matrix */
    947:         PetscCall(DMDAGetElementEqnums3D_up(u_eqn, p_eqn, ei, ej, ek));

    949:         for (n = 0; n < NODES_PER_EL; n++) {
    950:           if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) {
    951:             _ZERO_ROWCOL_i(Ae, 3 * n);
    952:             _ZERO_ROW_i(Ge, 3 * n);
    953:             _ZERO_COL_i(De, 3 * n);
    954:           }

    956:           if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) {
    957:             _ZERO_ROWCOL_i(Ae, 3 * n + 1);
    958:             _ZERO_ROW_i(Ge, 3 * n + 1);
    959:             _ZERO_COL_i(De, 3 * n + 1);
    960:           }

    962:           if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) {
    963:             _ZERO_ROWCOL_i(Ae, 3 * n + 2);
    964:             _ZERO_ROW_i(Ge, 3 * n + 2);
    965:             _ZERO_COL_i(De, 3 * n + 2);
    966:           }
    967:         }
    968:         PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
    969:         PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
    970:         PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
    971:         PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
    972:       }
    973:     }
    974:   }
    975:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    976:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    978:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    979:   PetscFunctionReturn(PETSC_SUCCESS);
    980: }

    982: static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, CellProperties cell_properties)
    983: {
    984:   DM                      cda;
    985:   Vec                     coords;
    986:   DMDACoor3d           ***_coords;
    987:   MatStencil              u_eqn[NODES_PER_EL * U_DOFS];
    988:   MatStencil              p_eqn[NODES_PER_EL * P_DOFS];
    989:   PetscInt                sex, sey, sez, mx, my, mz;
    990:   PetscInt                ei, ej, ek;
    991:   PetscScalar             Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
    992:   PetscScalar             Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
    993:   PetscScalar             De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
    994:   PetscScalar             Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
    995:   PetscScalar             el_coords[NODES_PER_EL * NSD];
    996:   GaussPointCoefficients *props;
    997:   PetscScalar            *prop_eta;
    998:   PetscInt                n, M, N, P;

    1000:   PetscFunctionBeginUser;
    1001:   PetscCall(DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    1002:   /* setup for coords */
    1003:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    1004:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    1005:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    1007:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
    1008:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
    1009:   for (ek = sez; ek < sez + mz; ek++) {
    1010:     for (ej = sey; ej < sey + my; ej++) {
    1011:       for (ei = sex; ei < sex + mx; ei++) {
    1012:         /* get coords for the element */
    1013:         PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
    1014:         /* get cell properties */
    1015:         PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &props));
    1016:         /* get coefficients for the element */
    1017:         prop_eta = props->eta;

    1019:         /* initialise element stiffness matrix */
    1020:         PetscCall(PetscMemzero(Ae, sizeof(Ae)));
    1021:         PetscCall(PetscMemzero(Ge, sizeof(Ge)));
    1022:         PetscCall(PetscMemzero(De, sizeof(De)));
    1023:         PetscCall(PetscMemzero(Ce, sizeof(Ce)));

    1025:         /* form element stiffness matrix */
    1026:         FormStressOperatorQ13D(Ae, el_coords, prop_eta);
    1027:         FormGradientOperatorQ13D(Ge, el_coords);
    1028:         /* FormDivergenceOperatorQ13D(De,el_coords); */
    1029:         FormScaledMassMatrixOperatorQ13D(Ce, el_coords, prop_eta);

    1031:         /* insert element matrix into global matrix */
    1032:         PetscCall(DMDAGetElementEqnums3D_up(u_eqn, p_eqn, ei, ej, ek));

    1034:         for (n = 0; n < NODES_PER_EL; n++) {
    1035:           if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) {
    1036:             _ZERO_ROWCOL_i(Ae, 3 * n);
    1037:             _ZERO_ROW_i(Ge, 3 * n);
    1038:             _ZERO_COL_i(De, 3 * n);
    1039:           }

    1041:           if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) {
    1042:             _ZERO_ROWCOL_i(Ae, 3 * n + 1);
    1043:             _ZERO_ROW_i(Ge, 3 * n + 1);
    1044:             _ZERO_COL_i(De, 3 * n + 1);
    1045:           }

    1047:           if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) {
    1048:             _ZERO_ROWCOL_i(Ae, 3 * n + 2);
    1049:             _ZERO_ROW_i(Ge, 3 * n + 2);
    1050:             _ZERO_COL_i(De, 3 * n + 2);
    1051:           }
    1052:         }
    1053:         PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
    1054:         PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
    1055:         /*PetscCall(MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES));*/
    1056:         PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
    1057:       }
    1058:     }
    1059:   }
    1060:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    1061:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    1063:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1064:   PetscFunctionReturn(PETSC_SUCCESS);
    1065: }

    1067: static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, CellProperties cell_properties)
    1068: {
    1069:   DM                      cda;
    1070:   Vec                     coords;
    1071:   DMDACoor3d           ***_coords;
    1072:   MatStencil              u_eqn[NODES_PER_EL * U_DOFS];
    1073:   MatStencil              p_eqn[NODES_PER_EL * P_DOFS];
    1074:   PetscInt                sex, sey, sez, mx, my, mz;
    1075:   PetscInt                ei, ej, ek;
    1076:   PetscScalar             Fe[NODES_PER_EL * U_DOFS];
    1077:   PetscScalar             He[NODES_PER_EL * P_DOFS];
    1078:   PetscScalar             el_coords[NODES_PER_EL * NSD];
    1079:   GaussPointCoefficients *props;
    1080:   PetscScalar            *prop_fx, *prop_fy, *prop_fz, *prop_hc;
    1081:   Vec                     local_F;
    1082:   StokesDOF            ***ff;
    1083:   PetscInt                n, M, N, P;

    1085:   PetscFunctionBeginUser;
    1086:   PetscCall(DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    1087:   /* setup for coords */
    1088:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    1089:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    1090:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    1092:   /* get access to the vector */
    1093:   PetscCall(DMGetLocalVector(stokes_da, &local_F));
    1094:   PetscCall(VecZeroEntries(local_F));
    1095:   PetscCall(DMDAVecGetArray(stokes_da, local_F, &ff));
    1096:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
    1097:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
    1098:   for (ek = sez; ek < sez + mz; ek++) {
    1099:     for (ej = sey; ej < sey + my; ej++) {
    1100:       for (ei = sex; ei < sex + mx; ei++) {
    1101:         /* get coords for the element */
    1102:         PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
    1103:         /* get cell properties */
    1104:         PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &props));
    1105:         /* get coefficients for the element */
    1106:         prop_fx = props->fx;
    1107:         prop_fy = props->fy;
    1108:         prop_fz = props->fz;
    1109:         prop_hc = props->hc;

    1111:         /* initialise element stiffness matrix */
    1112:         PetscCall(PetscMemzero(Fe, sizeof(Fe)));
    1113:         PetscCall(PetscMemzero(He, sizeof(He)));

    1115:         /* form element stiffness matrix */
    1116:         FormMomentumRhsQ13D(Fe, el_coords, prop_fx, prop_fy, prop_fz);
    1117:         FormContinuityRhsQ13D(He, el_coords, prop_hc);

    1119:         /* insert element matrix into global matrix */
    1120:         PetscCall(DMDAGetElementEqnums3D_up(u_eqn, p_eqn, ei, ej, ek));

    1122:         for (n = 0; n < NODES_PER_EL; n++) {
    1123:           if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) Fe[3 * n] = 0.0;

    1125:           if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) Fe[3 * n + 1] = 0.0;

    1127:           if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) Fe[3 * n + 2] = 0.0;
    1128:         }

    1130:         PetscCall(DMDASetValuesLocalStencil3D_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He));
    1131:       }
    1132:     }
    1133:   }
    1134:   PetscCall(DMDAVecRestoreArray(stokes_da, local_F, &ff));
    1135:   PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
    1136:   PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
    1137:   PetscCall(DMRestoreLocalVector(stokes_da, &local_F));

    1139:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1140:   PetscFunctionReturn(PETSC_SUCCESS);
    1141: }

    1143: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta, PetscReal *MX, PetscReal *MY, PetscReal *MZ)
    1144: {
    1145:   *theta = 0.0;
    1146:   *MX    = 2.0 * PETSC_PI;
    1147:   *MY    = 2.0 * PETSC_PI;
    1148:   *MZ    = 2.0 * PETSC_PI;
    1149: }
    1150: static void evaluate_MS_FrankKamentski(PetscReal pos[], PetscReal v[], PetscReal *p, PetscReal *eta, PetscReal Fm[], PetscReal *Fc)
    1151: {
    1152:   PetscReal x, y, z;
    1153:   PetscReal theta, MX, MY, MZ;

    1155:   evaluate_MS_FrankKamentski_constants(&theta, &MX, &MY, &MZ);
    1156:   x = pos[0];
    1157:   y = pos[1];
    1158:   z = pos[2];
    1159:   if (v) {
    1160:     /*
    1161:      v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
    1162:      v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
    1163:      v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
    1164:      */
    1165:     /*
    1166:      v[0] = PetscSinReal(PETSC_PI*x);
    1167:      v[1] = PetscSinReal(PETSC_PI*y);
    1168:      v[2] = PetscSinReal(PETSC_PI*z);
    1169:      */
    1170:     v[0] = PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x);
    1171:     v[1] = z * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y);
    1172:     v[2] = PetscPowRealInt(z, 2) * (PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 - PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2) - PETSC_PI * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) / 4;
    1173:   }
    1174:   if (p) *p = PetscPowRealInt(x, 2) + PetscPowRealInt(y, 2) + PetscPowRealInt(z, 2);
    1175:   if (eta) {
    1176:     /* eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
    1177:     *eta = 1.0;
    1178:   }
    1179:   if (Fm) {
    1180:     /*
    1181:      Fm[0] = -2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(PETSC_PI*x) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) - 0.2*PETSC_PI*MX*theta*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscCosReal(MZ*z)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 2.0*(3.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 3.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 1.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
    1182:      Fm[1] = -2*y - 0.2*MX*theta*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 2*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));
    1183:      Fm[2] = -2*z + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(2.0*PETSC_PI*z) - 0.2*MX*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 0.4*PETSC_PI*MZ*theta*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(MX*x)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-3.0*x*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(-3.0*y*PetscSinReal(2.0*PETSC_PI*z) - 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))  ;
    1184:      */
    1185:     /*
    1186:      Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
    1187:      Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
    1188:      Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
    1189:      */
    1190:     /*
    1191:      Fm[0] = -2*x + PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 6.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 6.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z) - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 2.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x) ;
    1192:      Fm[1] = -2*y - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z) + 2.0*z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 6.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 4.0*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
    1193:      Fm[2] = -2*z - 6.0*x*PetscSinReal(2.0*PETSC_PI*z) - 6.0*y*PetscSinReal(2.0*PETSC_PI*z) - PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z) + 3.0*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;
    1194:      */
    1195:     Fm[0] = -2 * x + 2 * z * (PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(PETSC_PI * y) * PetscExpReal(-y) * PetscSinReal(2.0 * PETSC_PI * x) - 1.0 * PETSC_PI * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) * PetscSinReal(2.0 * PETSC_PI * x)) + PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) + 6.0 * z * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) - 1.0 * PetscPowRealInt(PETSC_PI, 2) * PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) + 2.0 * PETSC_PI * z * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) * PetscSinReal(2.0 * PETSC_PI * x) - 2.0 * z * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(PETSC_PI * y) * PetscExpReal(-y) * PetscSinReal(2.0 * PETSC_PI * x);
    1196:     Fm[1] = -2 * y + 2 * z * (-PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y)) + 2.0 * z * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) - 6.0 * z * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) - 4.0 * PETSC_PI * z * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y);
    1197:     Fm[2] = -2 * z + PetscPowRealInt(z, 2) * (-2.0 * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) + 2.0 * PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y)) + PetscPowRealInt(z, 2) * (PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 - 3 * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2 - 3 * PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2) + 1.0 * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) + 0.25 * PetscPowRealInt(PETSC_PI, 3) * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 0.25 * PETSC_PI * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 3.0 * PETSC_PI * PetscPowRealInt(z, 2) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 1.0 * PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y);
    1198:   }
    1199:   if (Fc) {
    1200:     /* Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;*/
    1201:     /* Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
    1202:     /* Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);*/
    1203:     *Fc = 0.0;
    1204:   }
    1205: }

    1207: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx, PetscInt my, PetscInt mz, DM *_da, Vec *_X)
    1208: {
    1209:   DM            da, cda;
    1210:   Vec           X;
    1211:   StokesDOF  ***_stokes;
    1212:   Vec           coords;
    1213:   DMDACoor3d ***_coords;
    1214:   PetscInt      si, sj, sk, ei, ej, ek, i, j, k;

    1216:   PetscFunctionBeginUser;
    1217:   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 4, 1, NULL, NULL, NULL, &da));
    1218:   PetscCall(DMSetFromOptions(da));
    1219:   PetscCall(DMSetUp(da));
    1220:   PetscCall(DMDASetFieldName(da, 0, "anlytic_Vx"));
    1221:   PetscCall(DMDASetFieldName(da, 1, "anlytic_Vy"));
    1222:   PetscCall(DMDASetFieldName(da, 2, "anlytic_Vz"));
    1223:   PetscCall(DMDASetFieldName(da, 3, "analytic_P"));

    1225:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));

    1227:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1228:   PetscCall(DMGetCoordinateDM(da, &cda));
    1229:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    1231:   PetscCall(DMCreateGlobalVector(da, &X));
    1232:   PetscCall(DMDAVecGetArray(da, X, &_stokes));

    1234:   PetscCall(DMDAGetCorners(da, &si, &sj, &sk, &ei, &ej, &ek));
    1235:   for (k = sk; k < sk + ek; k++) {
    1236:     for (j = sj; j < sj + ej; j++) {
    1237:       for (i = si; i < si + ei; i++) {
    1238:         PetscReal pos[NSD], pressure, vel[NSD];

    1240:         pos[0] = PetscRealPart(_coords[k][j][i].x);
    1241:         pos[1] = PetscRealPart(_coords[k][j][i].y);
    1242:         pos[2] = PetscRealPart(_coords[k][j][i].z);

    1244:         evaluate_MS_FrankKamentski(pos, vel, &pressure, NULL, NULL, NULL);

    1246:         _stokes[k][j][i].u_dof = vel[0];
    1247:         _stokes[k][j][i].v_dof = vel[1];
    1248:         _stokes[k][j][i].w_dof = vel[2];
    1249:         _stokes[k][j][i].p_dof = pressure;
    1250:       }
    1251:     }
    1252:   }
    1253:   PetscCall(DMDAVecRestoreArray(da, X, &_stokes));
    1254:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

    1256:   *_da = da;
    1257:   *_X  = X;
    1258:   PetscFunctionReturn(PETSC_SUCCESS);
    1259: }

    1261: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da, Vec X, Vec X_analytic)
    1262: {
    1263:   DM            cda;
    1264:   Vec           coords, X_analytic_local, X_local;
    1265:   DMDACoor3d ***_coords;
    1266:   PetscInt      sex, sey, sez, mx, my, mz;
    1267:   PetscInt      ei, ej, ek;
    1268:   PetscScalar   el_coords[NODES_PER_EL * NSD];
    1269:   StokesDOF  ***stokes_analytic, ***stokes;
    1270:   StokesDOF     stokes_analytic_e[NODES_PER_EL], stokes_e[NODES_PER_EL];

    1272:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    1273:   PetscScalar Ni_p[NODES_PER_EL];
    1274:   PetscInt    ngp;
    1275:   PetscScalar gp_xi[GAUSS_POINTS][NSD];
    1276:   PetscScalar gp_weight[GAUSS_POINTS];
    1277:   PetscInt    p, i;
    1278:   PetscScalar J_p, fac;
    1279:   PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
    1280:   PetscScalar tint_p_ms, tint_p, int_p_ms, int_p;
    1281:   PetscInt    M;
    1282:   PetscReal   xymin[NSD], xymax[NSD];

    1284:   PetscFunctionBeginUser;
    1285:   /* define quadrature rule */
    1286:   ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    1288:   /* setup for coords */
    1289:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    1290:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    1291:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    1293:   /* setup for analytic */
    1294:   PetscCall(DMCreateLocalVector(stokes_da, &X_analytic_local));
    1295:   PetscCall(DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
    1296:   PetscCall(DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
    1297:   PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));

    1299:   /* setup for solution */
    1300:   PetscCall(DMCreateLocalVector(stokes_da, &X_local));
    1301:   PetscCall(DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local));
    1302:   PetscCall(DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local));
    1303:   PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));

    1305:   PetscCall(DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    1306:   PetscCall(DMGetBoundingBox(stokes_da, xymin, xymax));

    1308:   h = (xymax[0] - xymin[0]) / ((PetscReal)(M - 1));

    1310:   tp_L2 = tu_L2 = tu_H1 = 0.0;
    1311:   tint_p_ms = tint_p = 0.0;

    1313:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
    1314:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
    1315:   for (ek = sez; ek < sez + mz; ek++) {
    1316:     for (ej = sey; ej < sey + my; ej++) {
    1317:       for (ei = sex; ei < sex + mx; ei++) {
    1318:         /* get coords for the element */
    1319:         PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
    1320:         PetscCall(StokesDAGetNodalFields3D(stokes, ei, ej, ek, stokes_e));
    1321:         PetscCall(StokesDAGetNodalFields3D(stokes_analytic, ei, ej, ek, stokes_analytic_e));

    1323:         /* evaluate integral */
    1324:         for (p = 0; p < ngp; p++) {
    1325:           ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
    1326:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    1327:           ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, el_coords, &J_p);
    1328:           fac = gp_weight[p] * J_p;

    1330:           for (i = 0; i < NODES_PER_EL; i++) {
    1331:             tint_p_ms = tint_p_ms + fac * Ni_p[i] * stokes_analytic_e[i].p_dof;
    1332:             tint_p    = tint_p + fac * Ni_p[i] * stokes_e[i].p_dof;
    1333:           }
    1334:         }
    1335:       }
    1336:     }
    1337:   }
    1338:   PetscCallMPI(MPIU_Allreduce(&tint_p_ms, &int_p_ms, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    1339:   PetscCallMPI(MPIU_Allreduce(&tint_p, &int_p, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));

    1341:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\\int P dv %1.4e (h)  %1.4e (ms)\n", (double)PetscRealPart(int_p), (double)PetscRealPart(int_p_ms)));

    1343:   /* remove mine and add manufacture one */
    1344:   PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
    1345:   PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));

    1347:   {
    1348:     PetscInt     k, L, dof;
    1349:     PetscScalar *fields;
    1350:     PetscCall(DMDAGetInfo(stokes_da, 0, 0, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));

    1352:     PetscCall(VecGetLocalSize(X_local, &L));
    1353:     PetscCall(VecGetArray(X_local, &fields));
    1354:     for (k = 0; k < L / dof; k++) fields[dof * k + 3] = fields[dof * k + 3] - int_p + int_p_ms;
    1355:     PetscCall(VecRestoreArray(X_local, &fields));

    1357:     PetscCall(VecGetLocalSize(X, &L));
    1358:     PetscCall(VecGetArray(X, &fields));
    1359:     for (k = 0; k < L / dof; k++) fields[dof * k + 3] = fields[dof * k + 3] - int_p + int_p_ms;
    1360:     PetscCall(VecRestoreArray(X, &fields));
    1361:   }

    1363:   PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));
    1364:   PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));

    1366:   for (ek = sez; ek < sez + mz; ek++) {
    1367:     for (ej = sey; ej < sey + my; ej++) {
    1368:       for (ei = sex; ei < sex + mx; ei++) {
    1369:         /* get coords for the element */
    1370:         PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
    1371:         PetscCall(StokesDAGetNodalFields3D(stokes, ei, ej, ek, stokes_e));
    1372:         PetscCall(StokesDAGetNodalFields3D(stokes_analytic, ei, ej, ek, stokes_analytic_e));

    1374:         /* evaluate integral */
    1375:         p_e_L2 = 0.0;
    1376:         u_e_L2 = 0.0;
    1377:         u_e_H1 = 0.0;
    1378:         for (p = 0; p < ngp; p++) {
    1379:           ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
    1380:           ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
    1381:           ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, el_coords, &J_p);
    1382:           fac = gp_weight[p] * J_p;

    1384:           for (i = 0; i < NODES_PER_EL; i++) {
    1385:             PetscScalar u_error, v_error, w_error;

    1387:             p_e_L2 = p_e_L2 + fac * Ni_p[i] * (stokes_e[i].p_dof - stokes_analytic_e[i].p_dof) * (stokes_e[i].p_dof - stokes_analytic_e[i].p_dof);

    1389:             u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
    1390:             v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
    1391:             w_error = stokes_e[i].w_dof - stokes_analytic_e[i].w_dof;
    1392:             u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error + w_error * w_error);

    1394:             u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error                                                                                                   /* du/dx */
    1395:                                      + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error + GNx_p[2][i] * u_error * GNx_p[2][i] * u_error + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error /* dv/dx */
    1396:                                      + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error + GNx_p[2][i] * v_error * GNx_p[2][i] * v_error + GNx_p[0][i] * w_error * GNx_p[0][i] * w_error /* dw/dx */
    1397:                                      + GNx_p[1][i] * w_error * GNx_p[1][i] * w_error + GNx_p[2][i] * w_error * GNx_p[2][i] * w_error);
    1398:           }
    1399:         }

    1401:         tp_L2 += p_e_L2;
    1402:         tu_L2 += u_e_L2;
    1403:         tu_H1 += u_e_H1;
    1404:       }
    1405:     }
    1406:   }
    1407:   PetscCallMPI(MPIU_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    1408:   PetscCallMPI(MPIU_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    1409:   PetscCallMPI(MPIU_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    1410:   p_L2 = PetscSqrtScalar(p_L2);
    1411:   u_L2 = PetscSqrtScalar(u_L2);
    1412:   u_H1 = PetscSqrtScalar(u_H1);

    1414:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%1.4e   %1.4e   %1.4e   %1.4e \n", (double)PetscRealPart(h), (double)PetscRealPart(p_L2), (double)PetscRealPart(u_L2), (double)PetscRealPart(u_H1)));

    1416:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

    1418:   PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
    1419:   PetscCall(VecDestroy(&X_analytic_local));
    1420:   PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
    1421:   PetscCall(VecDestroy(&X_local));
    1422:   PetscFunctionReturn(PETSC_SUCCESS);
    1423: }

    1425: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da, Vec FIELD, const char file_prefix[])
    1426: {
    1427:   char         vtk_filename[PETSC_MAX_PATH_LEN];
    1428:   PetscMPIInt  rank;
    1429:   MPI_Comm     comm;
    1430:   FILE        *vtk_fp     = NULL;
    1431:   const char  *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
    1432:   PetscInt     si, sj, sk, nx, ny, nz, i;
    1433:   PetscInt     f, n_fields, N;
    1434:   DM           cda;
    1435:   Vec          coords;
    1436:   Vec          l_FIELD;
    1437:   PetscScalar *_L_FIELD;
    1438:   PetscInt     memory_offset;
    1439:   PetscScalar *buffer;

    1441:   PetscFunctionBeginUser;
    1442:   /* create file name */
    1443:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
    1444:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1445:   PetscCall(PetscSNPrintf(vtk_filename, sizeof(vtk_filename), "subdomain-%s-p%1.4d.vts", file_prefix, rank));

    1447:   /* open file and write header */
    1448:   vtk_fp = fopen(vtk_filename, "w");
    1449:   PetscCheck(vtk_fp, PETSC_COMM_SELF, PETSC_ERR_SYS, "Cannot open file = %s ", vtk_filename);

    1451:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<?xml version=\"1.0\"?>\n"));

    1453:   /* coords */
    1454:   PetscCall(DMDAGetGhostCorners(da, &si, &sj, &sk, &nx, &ny, &nz));
    1455:   N = nx * ny * nz;

    1457:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order));
    1458:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  <StructuredGrid WholeExtent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", si, si + nx - 1, sj, sj + ny - 1, sk, sk + nz - 1));
    1459:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", si, si + nx - 1, sj, sj + ny - 1, sk, sk + nz - 1));

    1461:   memory_offset = 0;

    1463:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      <CellData></CellData>\n"));

    1465:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      <Points>\n"));

    1467:   /* copy coordinates */
    1468:   PetscCall(DMGetCoordinateDM(da, &cda));
    1469:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1470:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "        <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n", memory_offset));
    1471:   memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar) * N * 3;

    1473:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      </Points>\n"));

    1475:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      <PointData Scalars=\" "));
    1476:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_fields, 0, 0, 0, 0, 0));
    1477:   for (f = 0; f < n_fields; f++) {
    1478:     const char *field_name;
    1479:     PetscCall(DMDAGetFieldName(da, f, &field_name));
    1480:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "%s ", field_name));
    1481:   }
    1482:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "\">\n"));

    1484:   for (f = 0; f < n_fields; f++) {
    1485:     const char *field_name;

    1487:     PetscCall(DMDAGetFieldName(da, f, &field_name));
    1488:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "        <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%" PetscInt_FMT "\"/>\n", field_name, memory_offset));
    1489:     memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar) * N;
    1490:   }

    1492:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      </PointData>\n"));
    1493:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    </Piece>\n"));
    1494:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  </StructuredGrid>\n"));

    1496:   PetscCall(PetscMalloc1(N, &buffer));
    1497:   PetscCall(DMGetLocalVector(da, &l_FIELD));
    1498:   PetscCall(DMGlobalToLocalBegin(da, FIELD, INSERT_VALUES, l_FIELD));
    1499:   PetscCall(DMGlobalToLocalEnd(da, FIELD, INSERT_VALUES, l_FIELD));
    1500:   PetscCall(VecGetArray(l_FIELD, &_L_FIELD));

    1502:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  <AppendedData encoding=\"raw\">\n"));
    1503:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "_"));

    1505:   /* write coordinates */
    1506:   {
    1507:     int          length = (int)(sizeof(PetscScalar) * N * 3);
    1508:     PetscScalar *allcoords;

    1510:     fwrite(&length, sizeof(int), 1, vtk_fp);
    1511:     PetscCall(VecGetArray(coords, &allcoords));
    1512:     fwrite(allcoords, sizeof(PetscScalar), 3 * N, vtk_fp);
    1513:     PetscCall(VecRestoreArray(coords, &allcoords));
    1514:   }
    1515:   /* write fields */
    1516:   for (f = 0; f < n_fields; f++) {
    1517:     int length = (int)(sizeof(PetscScalar) * N);
    1518:     fwrite(&length, sizeof(int), 1, vtk_fp);
    1519:     /* load */
    1520:     for (i = 0; i < N; i++) buffer[i] = _L_FIELD[n_fields * i + f];

    1522:     /* write */
    1523:     fwrite(buffer, sizeof(PetscScalar), N, vtk_fp);
    1524:   }
    1525:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "\n  </AppendedData>\n"));

    1527:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "</VTKFile>\n"));

    1529:   PetscCall(PetscFree(buffer));
    1530:   PetscCall(VecRestoreArray(l_FIELD, &_L_FIELD));
    1531:   PetscCall(DMRestoreLocalVector(da, &l_FIELD));

    1533:   if (vtk_fp) {
    1534:     fclose(vtk_fp);
    1535:     vtk_fp = NULL;
    1536:   }
    1537:   PetscFunctionReturn(PETSC_SUCCESS);
    1538: }

    1540: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp, PetscInt indent_level, DM da, const char local_file_prefix[])
    1541: {
    1542:   PetscMPIInt     size, rank;
    1543:   MPI_Comm        comm;
    1544:   const PetscInt *lx, *ly, *lz;
    1545:   PetscInt        M, N, P, pM, pN, pP, sum, *olx, *oly, *olz;
    1546:   PetscInt       *osx, *osy, *osz, *oex, *oey, *oez;
    1547:   PetscInt        i, j, k, II, stencil;

    1549:   PetscFunctionBeginUser;
    1550:   /* create file name */
    1551:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
    1552:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1553:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    1555:   PetscCall(DMDAGetInfo(da, 0, &M, &N, &P, &pM, &pN, &pP, 0, &stencil, 0, 0, 0, 0));
    1556:   PetscCall(DMDAGetOwnershipRanges(da, &lx, &ly, &lz));

    1558:   /* generate start,end list */
    1559:   PetscCall(PetscMalloc1(pM + 1, &olx));
    1560:   PetscCall(PetscMalloc1(pN + 1, &oly));
    1561:   PetscCall(PetscMalloc1(pP + 1, &olz));
    1562:   sum = 0;
    1563:   for (i = 0; i < pM; i++) {
    1564:     olx[i] = sum;
    1565:     sum    = sum + lx[i];
    1566:   }
    1567:   olx[pM] = sum;
    1568:   sum     = 0;
    1569:   for (i = 0; i < pN; i++) {
    1570:     oly[i] = sum;
    1571:     sum    = sum + ly[i];
    1572:   }
    1573:   oly[pN] = sum;
    1574:   sum     = 0;
    1575:   for (i = 0; i < pP; i++) {
    1576:     olz[i] = sum;
    1577:     sum    = sum + lz[i];
    1578:   }
    1579:   olz[pP] = sum;

    1581:   PetscCall(PetscMalloc1(pM, &osx));
    1582:   PetscCall(PetscMalloc1(pN, &osy));
    1583:   PetscCall(PetscMalloc1(pP, &osz));
    1584:   PetscCall(PetscMalloc1(pM, &oex));
    1585:   PetscCall(PetscMalloc1(pN, &oey));
    1586:   PetscCall(PetscMalloc1(pP, &oez));
    1587:   for (i = 0; i < pM; i++) {
    1588:     osx[i] = olx[i] - stencil;
    1589:     oex[i] = olx[i] + lx[i] + stencil;
    1590:     if (osx[i] < 0) osx[i] = 0;
    1591:     if (oex[i] > M) oex[i] = M;
    1592:   }

    1594:   for (i = 0; i < pN; i++) {
    1595:     osy[i] = oly[i] - stencil;
    1596:     oey[i] = oly[i] + ly[i] + stencil;
    1597:     if (osy[i] < 0) osy[i] = 0;
    1598:     if (oey[i] > M) oey[i] = N;
    1599:   }
    1600:   for (i = 0; i < pP; i++) {
    1601:     osz[i] = olz[i] - stencil;
    1602:     oez[i] = olz[i] + lz[i] + stencil;
    1603:     if (osz[i] < 0) osz[i] = 0;
    1604:     if (oez[i] > P) oez[i] = P;
    1605:   }

    1607:   for (k = 0; k < pP; k++) {
    1608:     for (j = 0; j < pN; j++) {
    1609:       for (i = 0; i < pM; i++) {
    1610:         char     name[PETSC_MAX_PATH_LEN];
    1611:         PetscInt procid = i + j * pM + k * pM * pN; /* convert proc(i,j,k) to pid */
    1612:         PetscCall(PetscSNPrintf(name, sizeof(name), "subdomain-%s-p%1.4" PetscInt_FMT ".vts", local_file_prefix, procid));
    1613:         for (II = 0; II < indent_level; II++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  "));

    1615:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\"      Source=\"%s\"/>\n", osx[i], oex[i] - 1, osy[j], oey[j] - 1, osz[k], oez[k] - 1, name));
    1616:       }
    1617:     }
    1618:   }
    1619:   PetscCall(PetscFree(olx));
    1620:   PetscCall(PetscFree(oly));
    1621:   PetscCall(PetscFree(olz));
    1622:   PetscCall(PetscFree(osx));
    1623:   PetscCall(PetscFree(osy));
    1624:   PetscCall(PetscFree(osz));
    1625:   PetscCall(PetscFree(oex));
    1626:   PetscCall(PetscFree(oey));
    1627:   PetscCall(PetscFree(oez));
    1628:   PetscFunctionReturn(PETSC_SUCCESS);
    1629: }

    1631: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da, const char file_prefix[], const char local_file_prefix[])
    1632: {
    1633:   MPI_Comm    comm;
    1634:   PetscMPIInt size, rank;
    1635:   char        vtk_filename[PETSC_MAX_PATH_LEN];
    1636:   FILE       *vtk_fp     = NULL;
    1637:   const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
    1638:   PetscInt    M, N, P, si, sj, sk, nx, ny, nz;
    1639:   PetscInt    i, dofs;

    1641:   PetscFunctionBeginUser;
    1642:   /* only rank-0 generates this file */
    1643:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
    1644:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1645:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    1647:   if (rank != 0) PetscFunctionReturn(PETSC_SUCCESS);

    1649:   /* create file name */
    1650:   PetscCall(PetscSNPrintf(vtk_filename, sizeof(vtk_filename), "%s.pvts", file_prefix));
    1651:   vtk_fp = fopen(vtk_filename, "w");
    1652:   PetscCheck(vtk_fp, PETSC_COMM_SELF, PETSC_ERR_SYS, "Cannot open file = %s ", vtk_filename);

    1654:   /* (VTK) generate pvts header */
    1655:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<?xml version=\"1.0\"?>\n"));

    1657:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order));

    1659:   /* define size of the nodal mesh based on the cell DM */
    1660:   PetscCall(DMDAGetInfo(da, 0, &M, &N, &P, 0, 0, 0, &dofs, 0, 0, 0, 0, 0));
    1661:   PetscCall(DMDAGetGhostCorners(da, &si, &sj, &sk, &nx, &ny, &nz));
    1662:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  <PStructuredGrid GhostLevel=\"1\" WholeExtent=\"%d %" PetscInt_FMT " %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, M - 1, 0, N - 1, 0, P - 1)); /* note overlap = 1 for Q1 */

    1664:   /* DUMP THE CELL REFERENCES */
    1665:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    <PCellData>\n"));
    1666:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    </PCellData>\n"));

    1668:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    <PPoints>\n"));
    1669:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n", NSD));
    1670:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    </PPoints>\n"));

    1672:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    <PPointData>\n"));
    1673:   for (i = 0; i < dofs; i++) {
    1674:     const char *fieldname;
    1675:     PetscCall(DMDAGetFieldName(da, i, &fieldname));
    1676:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "      <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n", fieldname));
    1677:   }
    1678:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "    </PPointData>\n"));

    1680:   /* write out the parallel information */
    1681:   PetscCall(DAViewVTK_write_PieceExtend(vtk_fp, 2, da, local_file_prefix));

    1683:   /* close the file */
    1684:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "  </PStructuredGrid>\n"));
    1685:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "</VTKFile>\n"));

    1687:   if (vtk_fp) {
    1688:     fclose(vtk_fp);
    1689:     vtk_fp = NULL;
    1690:   }
    1691:   PetscFunctionReturn(PETSC_SUCCESS);
    1692: }

    1694: PetscErrorCode DAView3DPVTS(DM da, Vec x, const char NAME[])
    1695: {
    1696:   char vts_filename[PETSC_MAX_PATH_LEN];
    1697:   char pvts_filename[PETSC_MAX_PATH_LEN];

    1699:   PetscFunctionBeginUser;
    1700:   PetscCall(PetscSNPrintf(vts_filename, sizeof(vts_filename), "%s-mesh", NAME));
    1701:   PetscCall(DAView_3DVTK_StructuredGrid_appended(da, x, vts_filename));

    1703:   PetscCall(PetscSNPrintf(pvts_filename, sizeof(pvts_filename), "%s-mesh", NAME));
    1704:   PetscCall(DAView_3DVTK_PStructuredGrid(da, pvts_filename, vts_filename));
    1705:   PetscFunctionReturn(PETSC_SUCCESS);
    1706: }

    1708: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp, PetscInt n, PetscReal rnorm, void *dummy)
    1709: {
    1710:   PetscReal norms[4];
    1711:   Vec       Br, v, w;
    1712:   Mat       A;

    1714:   PetscFunctionBeginUser;
    1715:   PetscCall(KSPGetOperators(ksp, &A, NULL));
    1716:   PetscCall(MatCreateVecs(A, &w, &v));

    1718:   PetscCall(KSPBuildResidual(ksp, v, w, &Br));

    1720:   PetscCall(VecStrideNorm(Br, 0, NORM_2, &norms[0]));
    1721:   PetscCall(VecStrideNorm(Br, 1, NORM_2, &norms[1]));
    1722:   PetscCall(VecStrideNorm(Br, 2, NORM_2, &norms[2]));
    1723:   PetscCall(VecStrideNorm(Br, 3, NORM_2, &norms[3]));

    1725:   PetscCall(VecDestroy(&v));
    1726:   PetscCall(VecDestroy(&w));

    1728:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%3" PetscInt_FMT " KSP Component U,V,W,P residual norm [ %1.12e, %1.12e, %1.12e, %1.12e ]\n", n, (double)norms[0], (double)norms[1], (double)norms[2], (double)norms[3]));
    1729:   PetscFunctionReturn(PETSC_SUCCESS);
    1730: }

    1732: static PetscErrorCode PCMGSetupViaCoarsen(PC pc, DM da_fine)
    1733: {
    1734:   PetscInt              nlevels, k;
    1735:   PETSC_UNUSED PetscInt finest;
    1736:   DM                   *da_list, *daclist;
    1737:   Mat                   R;

    1739:   PetscFunctionBeginUser;
    1740:   nlevels = 1;
    1741:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-levels", &nlevels, 0));

    1743:   PetscCall(PetscMalloc1(nlevels, &da_list));
    1744:   for (k = 0; k < nlevels; k++) da_list[k] = NULL;
    1745:   PetscCall(PetscMalloc1(nlevels, &daclist));
    1746:   for (k = 0; k < nlevels; k++) daclist[k] = NULL;

    1748:   /* finest grid is nlevels - 1 */
    1749:   finest     = nlevels - 1;
    1750:   daclist[0] = da_fine;
    1751:   PetscCall(PetscObjectReference((PetscObject)da_fine));
    1752:   PetscCall(DMCoarsenHierarchy(da_fine, nlevels - 1, &daclist[1]));
    1753:   for (k = 0; k < nlevels; k++) {
    1754:     da_list[k] = daclist[nlevels - 1 - k];
    1755:     PetscCall(DMDASetUniformCoordinates(da_list[k], 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
    1756:   }

    1758:   PetscCall(PCMGSetLevels(pc, nlevels, NULL));
    1759:   PetscCall(PCMGSetType(pc, PC_MG_MULTIPLICATIVE));
    1760:   PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_PMAT));

    1762:   for (k = 1; k < nlevels; k++) {
    1763:     PetscCall(DMCreateInterpolation(da_list[k - 1], da_list[k], &R, NULL));
    1764:     PetscCall(PCMGSetInterpolation(pc, k, R));
    1765:     PetscCall(MatDestroy(&R));
    1766:   }

    1768:   /* tidy up */
    1769:   for (k = 0; k < nlevels; k++) PetscCall(DMDestroy(&da_list[k]));
    1770:   PetscCall(PetscFree(da_list));
    1771:   PetscCall(PetscFree(daclist));
    1772:   PetscFunctionReturn(PETSC_SUCCESS);
    1773: }

    1775: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx, PetscInt my, PetscInt mz)
    1776: {
    1777:   DM             da_Stokes;
    1778:   PetscInt       u_dof, p_dof, dof, stencil_width;
    1779:   Mat            A, B;
    1780:   DM             vel_cda;
    1781:   Vec            vel_coords;
    1782:   PetscInt       p;
    1783:   Vec            f, X, X1;
    1784:   DMDACoor3d  ***_vel_coords;
    1785:   PetscInt       its;
    1786:   KSP            ksp_S;
    1787:   PetscInt       model_definition = 0;
    1788:   PetscInt       ei, ej, ek, sex, sey, sez, Imx, Jmy, Kmz;
    1789:   CellProperties cell_properties;
    1790:   PetscBool      write_output = PETSC_FALSE, resolve = PETSC_FALSE;

    1792:   PetscFunctionBeginUser;
    1793:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-resolve", &resolve, NULL));
    1794:   /* Generate the da for velocity and pressure */
    1795:   /* Num nodes in each direction is mx+1, my+1, mz+1 */
    1796:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
    1797:   p_dof         = P_DOFS; /* p - pressure */
    1798:   dof           = u_dof + p_dof;
    1799:   stencil_width = 1;
    1800:   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, NULL, &da_Stokes));
    1801:   PetscCall(DMSetMatType(da_Stokes, MATAIJ));
    1802:   PetscCall(DMSetFromOptions(da_Stokes));
    1803:   PetscCall(DMSetUp(da_Stokes));
    1804:   PetscCall(DMDASetFieldName(da_Stokes, 0, "Vx"));
    1805:   PetscCall(DMDASetFieldName(da_Stokes, 1, "Vy"));
    1806:   PetscCall(DMDASetFieldName(da_Stokes, 2, "Vz"));
    1807:   PetscCall(DMDASetFieldName(da_Stokes, 3, "P"));

    1809:   /* unit box [0,1] x [0,1] x [0,1] */
    1810:   PetscCall(DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));

    1812:   /* create quadrature point info for PDE coefficients */
    1813:   PetscCall(CellPropertiesCreate(da_Stokes, &cell_properties));

    1815:   /* interpolate the coordinates to quadrature points */
    1816:   PetscCall(DMGetCoordinateDM(da_Stokes, &vel_cda));
    1817:   PetscCall(DMGetCoordinatesLocal(da_Stokes, &vel_coords));
    1818:   PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
    1819:   PetscCall(DMDAGetElementsCorners(da_Stokes, &sex, &sey, &sez));
    1820:   PetscCall(DMDAGetElementsSizes(da_Stokes, &Imx, &Jmy, &Kmz));
    1821:   for (ek = sez; ek < sez + Kmz; ek++) {
    1822:     for (ej = sey; ej < sey + Jmy; ej++) {
    1823:       for (ei = sex; ei < sex + Imx; ei++) {
    1824:         /* get coords for the element */
    1825:         PetscInt                ngp;
    1826:         PetscScalar             gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
    1827:         PetscScalar             el_coords[NSD * NODES_PER_EL];
    1828:         GaussPointCoefficients *cell;

    1830:         PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
    1831:         PetscCall(GetElementCoords3D(_vel_coords, ei, ej, ek, el_coords));
    1832:         ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);

    1834:         for (p = 0; p < GAUSS_POINTS; p++) {
    1835:           PetscScalar xi_p[NSD], Ni_p[NODES_PER_EL];
    1836:           PetscScalar gp_x, gp_y, gp_z;
    1837:           PetscInt    n;

    1839:           xi_p[0] = gp_xi[p][0];
    1840:           xi_p[1] = gp_xi[p][1];
    1841:           xi_p[2] = gp_xi[p][2];
    1842:           ShapeFunctionQ13D_Evaluate(xi_p, Ni_p);

    1844:           gp_x = gp_y = gp_z = 0.0;
    1845:           for (n = 0; n < NODES_PER_EL; n++) {
    1846:             gp_x = gp_x + Ni_p[n] * el_coords[NSD * n];
    1847:             gp_y = gp_y + Ni_p[n] * el_coords[NSD * n + 1];
    1848:             gp_z = gp_z + Ni_p[n] * el_coords[NSD * n + 2];
    1849:           }
    1850:           cell->gp_coords[NSD * p]     = gp_x;
    1851:           cell->gp_coords[NSD * p + 1] = gp_y;
    1852:           cell->gp_coords[NSD * p + 2] = gp_z;
    1853:         }
    1854:       }
    1855:     }
    1856:   }

    1858:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-model", &model_definition, NULL));

    1860:   switch (model_definition) {
    1861:   case 0: /* isoviscous */
    1862:     for (ek = sez; ek < sez + Kmz; ek++) {
    1863:       for (ej = sey; ej < sey + Jmy; ej++) {
    1864:         for (ei = sex; ei < sex + Imx; ei++) {
    1865:           GaussPointCoefficients *cell;

    1867:           PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
    1868:           for (p = 0; p < GAUSS_POINTS; p++) {
    1869:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
    1870:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
    1871:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);

    1873:             cell->eta[p] = 1.0;

    1875:             cell->fx[p] = 0.0 * coord_x;
    1876:             cell->fy[p] = -PetscSinReal(2.2 * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
    1877:             cell->fz[p] = 0.0 * coord_z;
    1878:             cell->hc[p] = 0.0;
    1879:           }
    1880:         }
    1881:       }
    1882:     }
    1883:     break;

    1885:   case 1: /* manufactured */
    1886:     for (ek = sez; ek < sez + Kmz; ek++) {
    1887:       for (ej = sey; ej < sey + Jmy; ej++) {
    1888:         for (ei = sex; ei < sex + Imx; ei++) {
    1889:           PetscReal               eta, Fm[NSD], Fc, pos[NSD];
    1890:           GaussPointCoefficients *cell;

    1892:           PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
    1893:           for (p = 0; p < GAUSS_POINTS; p++) {
    1894:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
    1895:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
    1896:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);

    1898:             pos[0] = coord_x;
    1899:             pos[1] = coord_y;
    1900:             pos[2] = coord_z;

    1902:             evaluate_MS_FrankKamentski(pos, NULL, NULL, &eta, Fm, &Fc);
    1903:             cell->eta[p] = eta;

    1905:             cell->fx[p] = Fm[0];
    1906:             cell->fy[p] = Fm[1];
    1907:             cell->fz[p] = Fm[2];
    1908:             cell->hc[p] = Fc;
    1909:           }
    1910:         }
    1911:       }
    1912:     }
    1913:     break;

    1915:   case 2: /* solcx */
    1916:     for (ek = sez; ek < sez + Kmz; ek++) {
    1917:       for (ej = sey; ej < sey + Jmy; ej++) {
    1918:         for (ei = sex; ei < sex + Imx; ei++) {
    1919:           GaussPointCoefficients *cell;

    1921:           PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
    1922:           for (p = 0; p < GAUSS_POINTS; p++) {
    1923:             PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
    1924:             PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
    1925:             PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);

    1927:             cell->eta[p] = 1.0;

    1929:             cell->fx[p] = 0.0;
    1930:             cell->fy[p] = -PetscSinReal(3.0 * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
    1931:             cell->fz[p] = 0.0 * coord_z;
    1932:             cell->hc[p] = 0.0;
    1933:           }
    1934:         }
    1935:       }
    1936:     }
    1937:     break;

    1939:   case 3: /* sinker */
    1940:     for (ek = sez; ek < sez + Kmz; ek++) {
    1941:       for (ej = sey; ej < sey + Jmy; ej++) {
    1942:         for (ei = sex; ei < sex + Imx; ei++) {
    1943:           GaussPointCoefficients *cell;

    1945:           PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
    1946:           for (p = 0; p < GAUSS_POINTS; p++) {
    1947:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD * p]);
    1948:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD * p + 1]);
    1949:             PetscReal zp = PetscRealPart(cell->gp_coords[NSD * p + 2]);

    1951:             cell->eta[p] = 1.0e-2;
    1952:             cell->fx[p]  = 0.0;
    1953:             cell->fy[p]  = 0.0;
    1954:             cell->fz[p]  = 0.0;
    1955:             cell->hc[p]  = 0.0;

    1957:             if ((PetscAbs(xp - 0.5) < 0.2) && (PetscAbs(yp - 0.5) < 0.2) && (PetscAbs(zp - 0.5) < 0.2)) {
    1958:               cell->eta[p] = 1.0;
    1959:               cell->fz[p]  = 1.0;
    1960:             }
    1961:           }
    1962:         }
    1963:       }
    1964:     }
    1965:     break;

    1967:   case 4: /* subdomain jumps */
    1968:     for (ek = sez; ek < sez + Kmz; ek++) {
    1969:       for (ej = sey; ej < sey + Jmy; ej++) {
    1970:         for (ei = sex; ei < sex + Imx; ei++) {
    1971:           PetscReal               opts_mag, opts_eta0;
    1972:           PetscInt                px, py, pz;
    1973:           PetscBool               jump;
    1974:           PetscMPIInt             rr;
    1975:           GaussPointCoefficients *cell;

    1977:           opts_mag  = 1.0;
    1978:           opts_eta0 = 1.e-2;

    1980:           PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL));
    1981:           PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL));
    1982:           PetscCall(DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, &pz, NULL, NULL, NULL, NULL, NULL, NULL));
    1983:           rr = PetscGlobalRank % (px * py);
    1984:           if (px % 2) jump = (PetscBool)(rr % 2);
    1985:           else jump = (PetscBool)((rr / px) % 2 ? rr % 2 : !(rr % 2));
    1986:           rr = PetscGlobalRank / (px * py);
    1987:           if (rr % 2) jump = (PetscBool)!jump;
    1988:           PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
    1989:           for (p = 0; p < GAUSS_POINTS; p++) {
    1990:             PetscReal xp = PetscRealPart(cell->gp_coords[NSD * p]);
    1991:             PetscReal yp = PetscRealPart(cell->gp_coords[NSD * p + 1]);

    1993:             cell->eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
    1994:             cell->fx[p]  = 0.0;
    1995:             cell->fy[p]  = -PetscSinReal(2.2 * PETSC_PI * yp) * PetscCosReal(1.0 * PETSC_PI * xp);
    1996:             cell->fz[p]  = 0.0;
    1997:             cell->hc[p]  = 0.0;
    1998:           }
    1999:         }
    2000:       }
    2001:     }
    2002:     break;
    2003:   default:
    2004:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "No default model is supported. Choose either -model {0,1,2,3}");
    2005:   }

    2007:   PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));

    2009:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
    2010:   PetscCall(DMCreateMatrix(da_Stokes, &A));
    2011:   PetscCall(DMCreateMatrix(da_Stokes, &B));
    2012:   PetscCall(DMCreateGlobalVector(da_Stokes, &X));
    2013:   PetscCall(DMCreateGlobalVector(da_Stokes, &f));

    2015:   /* assemble A11 */
    2016:   PetscCall(MatZeroEntries(A));
    2017:   PetscCall(MatZeroEntries(B));
    2018:   PetscCall(VecZeroEntries(f));

    2020:   PetscCall(AssembleA_Stokes(A, da_Stokes, cell_properties));
    2021:   PetscCall(MatViewFromOptions(A, NULL, "-amat_view"));
    2022:   PetscCall(AssembleA_PCStokes(B, da_Stokes, cell_properties));
    2023:   PetscCall(MatViewFromOptions(B, NULL, "-bmat_view"));
    2024:   /* build force vector */
    2025:   PetscCall(AssembleF_Stokes(f, da_Stokes, cell_properties));

    2027:   /* SOLVE */
    2028:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_S));
    2029:   PetscCall(KSPSetOptionsPrefix(ksp_S, "stokes_"));
    2030:   PetscCall(KSPSetOperators(ksp_S, A, B));
    2031:   PetscCall(KSPSetFromOptions(ksp_S));

    2033:   {
    2034:     PC             pc;
    2035:     const PetscInt ufields[] = {0, 1, 2}, pfields[] = {3};
    2036:     PetscCall(KSPGetPC(ksp_S, &pc));
    2037:     PetscCall(PCFieldSplitSetBlockSize(pc, 4));
    2038:     PetscCall(PCFieldSplitSetFields(pc, "u", 3, ufields, ufields));
    2039:     PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
    2040:   }

    2042:   {
    2043:     PC        pc;
    2044:     PetscBool same = PETSC_FALSE;
    2045:     PetscCall(KSPGetPC(ksp_S, &pc));
    2046:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCMG, &same));
    2047:     if (same) PetscCall(PCMGSetupViaCoarsen(pc, da_Stokes));
    2048:   }

    2050:   {
    2051:     PC        pc;
    2052:     PetscBool same = PETSC_FALSE;
    2053:     PetscCall(KSPGetPC(ksp_S, &pc));
    2054:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same));
    2055:     if (same) PetscCall(KSPSetOperators(ksp_S, A, A));
    2056:   }

    2058:   {
    2059:     PetscBool stokes_monitor = PETSC_FALSE;
    2060:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-stokes_ksp_monitor_blocks", &stokes_monitor, 0));
    2061:     if (stokes_monitor) PetscCall(KSPMonitorSet(ksp_S, KSPMonitorStokesBlocks, NULL, NULL));
    2062:   }

    2064:   if (resolve) {
    2065:     /* Test changing matrix data structure and resolve */
    2066:     PetscCall(VecDuplicate(X, &X1));
    2067:   }

    2069:   PetscCall(KSPSolve(ksp_S, f, X));
    2070:   if (resolve) {
    2071:     Mat C;
    2072:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &C));
    2073:     PetscCall(KSPSetOperators(ksp_S, C, C));
    2074:     PetscCall(KSPSolve(ksp_S, f, X1));
    2075:     PetscCall(MatDestroy(&C));
    2076:     PetscCall(VecDestroy(&X1));
    2077:   }

    2079:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-write_pvts", &write_output, NULL));
    2080:   if (write_output) PetscCall(DAView3DPVTS(da_Stokes, X, "up"));
    2081:   {
    2082:     PetscBool flg = PETSC_FALSE;
    2083:     char      filename[PETSC_MAX_PATH_LEN];
    2084:     PetscCall(PetscOptionsGetString(NULL, NULL, "-write_binary", filename, sizeof(filename), &flg));
    2085:     if (flg) {
    2086:       PetscViewer viewer;
    2087:       /* PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer)); */
    2088:       PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, "ex42.vts", FILE_MODE_WRITE, &viewer));
    2089:       PetscCall(VecView(X, viewer));
    2090:       PetscCall(PetscViewerDestroy(&viewer));
    2091:     }
    2092:   }
    2093:   PetscCall(KSPGetIterationNumber(ksp_S, &its));

    2095:   /* verify */
    2096:   if (model_definition == 1) {
    2097:     DM  da_Stokes_analytic;
    2098:     Vec X_analytic;

    2100:     PetscCall(DMDACreateManufacturedSolution(mx, my, mz, &da_Stokes_analytic, &X_analytic));
    2101:     if (write_output) PetscCall(DAView3DPVTS(da_Stokes_analytic, X_analytic, "ms"));
    2102:     PetscCall(DMDAIntegrateErrors3D(da_Stokes_analytic, X, X_analytic));
    2103:     if (write_output) PetscCall(DAView3DPVTS(da_Stokes, X, "up2"));
    2104:     PetscCall(DMDestroy(&da_Stokes_analytic));
    2105:     PetscCall(VecDestroy(&X_analytic));
    2106:   }

    2108:   PetscCall(KSPDestroy(&ksp_S));
    2109:   PetscCall(VecDestroy(&X));
    2110:   PetscCall(VecDestroy(&f));
    2111:   PetscCall(MatDestroy(&A));
    2112:   PetscCall(MatDestroy(&B));

    2114:   PetscCall(CellPropertiesDestroy(&cell_properties));
    2115:   PetscCall(DMDestroy(&da_Stokes));
    2116:   PetscFunctionReturn(PETSC_SUCCESS);
    2117: }

    2119: int main(int argc, char **args)
    2120: {
    2121:   PetscInt mx, my, mz;

    2123:   PetscFunctionBeginUser;
    2124:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
    2125:   mx = my = mz = 10;
    2126:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
    2127:   my = mx;
    2128:   mz = mx;
    2129:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
    2130:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, NULL));
    2131:   PetscCall(solve_stokes_3d_coupled(mx, my, mz));
    2132:   PetscCall(PetscFinalize());
    2133:   return 0;
    2134: }

    2136: /*TEST

    2138:    test:
    2139:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type lu

    2141:    test:
    2142:       suffix: 2
    2143:       nsize: 3
    2144:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type redundant -stokes_redundant_pc_type lu

    2146:    test:
    2147:       suffix: bddc_stokes
    2148:       nsize: 6
    2149:       args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -stokes_pc_bddc_coarse_redundant_pc_type svd

    2151:    test:
    2152:       suffix: bddc_stokes_deluxe
    2153:       nsize: 6
    2154:       args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc

    2156:    test:
    2157:       requires: !single
    2158:       suffix: bddc_stokes_subdomainjump_deluxe
    2159:       nsize: 9
    2160:       args: -model 4 -jump_magnitude 4 -mx 6 -my 6 -mz 2 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc -stokes_pc_bddc_schur_layers 1

    2162:    test:
    2163:       requires: !single
    2164:       suffix: 3
    2165:       args: -stokes_ksp_converged_reason -stokes_pc_type fieldsplit -resolve

    2167:    test:
    2168:       suffix: tut_1
    2169:       nsize: 4
    2170:       requires: !single
    2171:       args: -stokes_ksp_monitor

    2173:    test:
    2174:       suffix: tut_2
    2175:       nsize: 4
    2176:       requires: !single
    2177:       args: -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur

    2179:    test:
    2180:       suffix: tut_3
    2181:       nsize: 4
    2182:       requires: !single
    2183:       args: -mx 20 -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur

    2185: TEST*/