Actual source code: ex43.c

      1: static char help[] = "Solves the incompressible, variable viscosity Stokes equation in 2d on the unit domain \n\
      2: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
      3: The models defined utilise free slip boundary conditions on all sides. \n\
      4: Options: \n"
      5:                      "\
      6:      -mx : Number of elements in the x-direction \n\
      7:      -my : Number of elements in the y-direction \n\
      8:      -o : Specify output filename for solution (will be petsc binary format or paraview format if the extension is .vts) \n\
      9:      -gnuplot : Output Gauss point coordinates, coefficients and u,p solution in gnuplot format \n\
     10:      -glvis : Visualizes coefficients and u,p solution through GLVIs (use -viewer_glvis_dmda_bs 2,1 to visualize velocity as a vector)\n\
     11:      -c_str : Indicates the structure of the coefficients to use \n"
     12:                      "\
     13:           -c_str 0 => Coefficient definition for an analytic solution with a vertical jump in viscosity at x = xc \n\
     14:                       This problem is driven by the forcing function f(x,y) = (0, sin(nz pi y)cos(pi x) \n\
     15:                          Parameters: \n\
     16:                               -solcx_eta0  : Viscosity to the left of the interface \n\
     17:                               -solcx_eta1  : Viscosity to the right of the interface \n\
     18:                               -solcx_xc    : Location of the interface \n\
     19:                               -solcx_nz    : Wavenumber in the y direction \n"
     20:                      "\
     21:           -c_str 1 => Coefficient definition for a dense rectangular blob located at the center of the domain \n\
     22:                          Parameters: \n\
     23:                               -sinker_eta0 : Viscosity of the background fluid \n\
     24:                               -sinker_eta1 : Viscosity of the blob \n\
     25:                               -sinker_dx   : Width of the blob \n\
     26:                               -sinker_dy   : Height of the blob \n"
     27:                      "\
     28:           -c_str 2 => Coefficient definition for a dense circular blob located at the center of the domain \n\
     29:                          Parameters: \n\
     30:                               -sinker_eta0 : Viscosity of the background fluid \n\
     31:                               -sinker_eta1 : Viscosity of the blob \n\
     32:                               -sinker_r    : Radius of the blob \n"
     33:                      "\
     34:           -c_str 3 => Coefficient definition for a dense circular and rectangular inclusion (located at the center of the domain) \n\
     35:                               -sinker_eta0 : Viscosity of the background fluid \n\
     36:                               -sinker_eta1 : Viscosity of the two inclusions \n\
     37:                               -sinker_r    : Radius of the circular inclusion \n\
     38:                               -sinker_c0x  : Origin (x-coord) of the circular inclusion \n\
     39:                               -sinker_c0y  : Origin (y-coord) of the circular inclusion \n\
     40:                               -sinker_dx   : Width of the rectangular inclusion \n\
     41:                               -sinker_dy   : Height of the rectangular inclusion \n\
     42:                               -sinker_phi  : Rotation angle of the rectangular inclusion \n"
     43:                      "\
     44:           -c_str 4 => Coefficient definition for checkerboard jumps aligned with the domain decomposition \n\
     45:                               -jump_eta0      : Viscosity for black subdomains \n\
     46:                               -jump_magnitude : Magnitude of jumps. White subdomains will have eta = eta0*10^magnitude \n\
     47:                               -jump_nz        : Wavenumber in the y direction for rhs \n"
     48:                      "\
     49:      -use_gp_coords : Evaluate the viscosity and force term at the global coordinates of each quadrature point \n\
     50:                       By default, the viscosity and force term are evaluated at the element center and applied as a constant over the entire element \n";

     52: /* Contributed by Dave May */

     54: #include <petscksp.h>
     55: #include <petscdm.h>
     56: #include <petscdmda.h>

     58: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
     59: #include "ex43-solcx.h"

     61: static PetscErrorCode DMDABCApplyFreeSlip(DM, Mat, Vec);

     63: #define NSD          2 /* number of spatial dimensions */
     64: #define NODES_PER_EL 4 /* nodes per element */
     65: #define U_DOFS       2 /* degrees of freedom per velocity node */
     66: #define P_DOFS       1 /* degrees of freedom per pressure node */
     67: #define GAUSS_POINTS 4

     69: /* Gauss point based evaluation 8+4+4+4 = 20 */
     70: typedef struct {
     71:   PetscScalar gp_coords[2 * GAUSS_POINTS];
     72:   PetscScalar eta[GAUSS_POINTS];
     73:   PetscScalar fx[GAUSS_POINTS];
     74:   PetscScalar fy[GAUSS_POINTS];
     75: } GaussPointCoefficients;

     77: typedef struct {
     78:   PetscScalar u_dof;
     79:   PetscScalar v_dof;
     80:   PetscScalar p_dof;
     81: } StokesDOF;

     83: static PetscErrorCode glvis_extract_eta(PetscObject oV, PetscInt nf, PetscObject oVf[], void *ctx)
     84: {
     85:   DM                       properties_da = (DM)ctx, stokes_da;
     86:   Vec                      V = (Vec)oV, *Vf = (Vec *)oVf;
     87:   GaussPointCoefficients **props;
     88:   PetscInt                 sex, sey, mx, my;
     89:   PetscInt                 ei, ej, p, cum;
     90:   PetscScalar             *array;

     92:   PetscFunctionBeginUser;
     93:   PetscCall(VecGetDM(Vf[0], &stokes_da));
     94:   PetscCall(DMDAVecGetArrayRead(properties_da, V, &props));
     95:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
     96:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
     97:   PetscCall(VecGetArray(Vf[0], &array));
     98:   cum = 0;
     99:   for (ej = sey; ej < sey + my; ej++) {
    100:     for (ei = sex; ei < sex + mx; ei++) {
    101:       for (p = 0; p < GAUSS_POINTS; p++) array[cum++] = props[ej][ei].eta[p];
    102:     }
    103:   }
    104:   PetscCall(VecRestoreArray(Vf[0], &array));
    105:   PetscCall(DMDAVecRestoreArrayRead(properties_da, V, &props));
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: /*
    110:  Element: Local basis function ordering
    111:  1-----2
    112:  |     |
    113:  |     |
    114:  0-----3
    115:  */
    116: static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
    117: {
    118:   PetscScalar xi  = _xi[0];
    119:   PetscScalar eta = _xi[1];

    121:   Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
    122:   Ni[1] = 0.25 * (1.0 - xi) * (1.0 + eta);
    123:   Ni[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
    124:   Ni[3] = 0.25 * (1.0 + xi) * (1.0 - eta);
    125: }

    127: static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
    128: {
    129:   PetscScalar xi  = _xi[0];
    130:   PetscScalar eta = _xi[1];

    132:   GNi[0][0] = -0.25 * (1.0 - eta);
    133:   GNi[0][1] = -0.25 * (1.0 + eta);
    134:   GNi[0][2] = 0.25 * (1.0 + eta);
    135:   GNi[0][3] = 0.25 * (1.0 - eta);

    137:   GNi[1][0] = -0.25 * (1.0 - xi);
    138:   GNi[1][1] = 0.25 * (1.0 - xi);
    139:   GNi[1][2] = 0.25 * (1.0 + xi);
    140:   GNi[1][3] = -0.25 * (1.0 + xi);
    141: }

    143: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
    144: {
    145:   PetscScalar J00, J01, J10, J11, J;
    146:   PetscScalar iJ00, iJ01, iJ10, iJ11;
    147:   PetscInt    i;

    149:   J00 = J01 = J10 = J11 = 0.0;
    150:   for (i = 0; i < NODES_PER_EL; i++) {
    151:     PetscScalar cx = coords[2 * i];
    152:     PetscScalar cy = coords[2 * i + 1];

    154:     J00 = J00 + GNi[0][i] * cx; /* J_xx = dx/dxi */
    155:     J01 = J01 + GNi[0][i] * cy; /* J_xy = dy/dxi */
    156:     J10 = J10 + GNi[1][i] * cx; /* J_yx = dx/deta */
    157:     J11 = J11 + GNi[1][i] * cy; /* J_yy = dy/deta */
    158:   }
    159:   J = (J00 * J11) - (J01 * J10);

    161:   iJ00 = J11 / J;
    162:   iJ01 = -J01 / J;
    163:   iJ10 = -J10 / J;
    164:   iJ11 = J00 / J;

    166:   for (i = 0; i < NODES_PER_EL; i++) {
    167:     GNx[0][i] = GNi[0][i] * iJ00 + GNi[1][i] * iJ01;
    168:     GNx[1][i] = GNi[0][i] * iJ10 + GNi[1][i] * iJ11;
    169:   }

    171:   if (det_J) *det_J = J;
    172: }

    174: static void ConstructGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
    175: {
    176:   *ngp         = 4;
    177:   gp_xi[0][0]  = -0.57735026919;
    178:   gp_xi[0][1]  = -0.57735026919;
    179:   gp_xi[1][0]  = -0.57735026919;
    180:   gp_xi[1][1]  = 0.57735026919;
    181:   gp_xi[2][0]  = 0.57735026919;
    182:   gp_xi[2][1]  = 0.57735026919;
    183:   gp_xi[3][0]  = 0.57735026919;
    184:   gp_xi[3][1]  = -0.57735026919;
    185:   gp_weight[0] = 1.0;
    186:   gp_weight[1] = 1.0;
    187:   gp_weight[2] = 1.0;
    188:   gp_weight[3] = 1.0;
    189: }

    191: /*
    192: i,j are the element indices
    193: The unknown is a vector quantity.
    194: The s[].c is used to indicate the degree of freedom.
    195: */
    196: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[], MatStencil s_p[], PetscInt i, PetscInt j)
    197: {
    198:   PetscFunctionBeginUser;
    199:   /* velocity */
    200:   /* node 0 */
    201:   s_u[0].i = i;
    202:   s_u[0].j = j;
    203:   s_u[0].c = 0; /* Vx0 */
    204:   s_u[1].i = i;
    205:   s_u[1].j = j;
    206:   s_u[1].c = 1; /* Vy0 */

    208:   /* node 1 */
    209:   s_u[2].i = i;
    210:   s_u[2].j = j + 1;
    211:   s_u[2].c = 0; /* Vx1 */
    212:   s_u[3].i = i;
    213:   s_u[3].j = j + 1;
    214:   s_u[3].c = 1; /* Vy1 */

    216:   /* node 2 */
    217:   s_u[4].i = i + 1;
    218:   s_u[4].j = j + 1;
    219:   s_u[4].c = 0; /* Vx2 */
    220:   s_u[5].i = i + 1;
    221:   s_u[5].j = j + 1;
    222:   s_u[5].c = 1; /* Vy2 */

    224:   /* node 3 */
    225:   s_u[6].i = i + 1;
    226:   s_u[6].j = j;
    227:   s_u[6].c = 0; /* Vx3 */
    228:   s_u[7].i = i + 1;
    229:   s_u[7].j = j;
    230:   s_u[7].c = 1; /* Vy3 */

    232:   /* pressure */
    233:   s_p[0].i = i;
    234:   s_p[0].j = j;
    235:   s_p[0].c = 2; /* P0 */
    236:   s_p[1].i = i;
    237:   s_p[1].j = j + 1;
    238:   s_p[1].c = 2; /* P1 */
    239:   s_p[2].i = i + 1;
    240:   s_p[2].j = j + 1;
    241:   s_p[2].c = 2; /* P2 */
    242:   s_p[3].i = i + 1;
    243:   s_p[3].j = j;
    244:   s_p[3].c = 2; /* P3 */
    245:   PetscFunctionReturn(PETSC_SUCCESS);
    246: }

    248: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da, PetscInt **_lx, PetscInt **_ly)
    249: {
    250:   PetscMPIInt  rank;
    251:   PetscInt     proc_I, proc_J;
    252:   PetscInt     cpu_x, cpu_y;
    253:   PetscInt     local_mx, local_my;
    254:   Vec          vlx, vly;
    255:   PetscInt    *LX, *LY, i;
    256:   PetscScalar *_a;
    257:   Vec          V_SEQ;
    258:   VecScatter   ctx;

    260:   PetscFunctionBeginUser;
    261:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

    263:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));

    265:   proc_J = rank / cpu_x;
    266:   proc_I = rank - cpu_x * proc_J;

    268:   PetscCall(PetscMalloc1(cpu_x, &LX));
    269:   PetscCall(PetscMalloc1(cpu_y, &LY));

    271:   PetscCall(DMDAGetElementsSizes(da, &local_mx, &local_my, NULL));
    272:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vlx));
    273:   PetscCall(VecSetSizes(vlx, PETSC_DECIDE, cpu_x));
    274:   PetscCall(VecSetFromOptions(vlx));

    276:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
    277:   PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
    278:   PetscCall(VecSetFromOptions(vly));

    280:   PetscCall(VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES));
    281:   PetscCall(VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES));
    282:   PetscCall(VecAssemblyBegin(vlx));
    283:   PetscCall(VecAssemblyEnd(vlx));
    284:   PetscCall(VecAssemblyBegin(vly));
    285:   PetscCall(VecAssemblyEnd(vly));

    287:   PetscCall(VecScatterCreateToAll(vlx, &ctx, &V_SEQ));
    288:   PetscCall(VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    289:   PetscCall(VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    290:   PetscCall(VecGetArray(V_SEQ, &_a));
    291:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
    292:   PetscCall(VecRestoreArray(V_SEQ, &_a));
    293:   PetscCall(VecScatterDestroy(&ctx));
    294:   PetscCall(VecDestroy(&V_SEQ));

    296:   PetscCall(VecScatterCreateToAll(vly, &ctx, &V_SEQ));
    297:   PetscCall(VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    298:   PetscCall(VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    299:   PetscCall(VecGetArray(V_SEQ, &_a));
    300:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
    301:   PetscCall(VecRestoreArray(V_SEQ, &_a));
    302:   PetscCall(VecScatterDestroy(&ctx));
    303:   PetscCall(VecDestroy(&V_SEQ));

    305:   *_lx = LX;
    306:   *_ly = LY;

    308:   PetscCall(VecDestroy(&vlx));
    309:   PetscCall(VecDestroy(&vly));
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

    313: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
    314: {
    315:   DM           cda;
    316:   Vec          coords;
    317:   DMDACoor2d **_coords;
    318:   PetscInt     si, sj, nx, ny, i, j;
    319:   FILE        *fp;
    320:   char         fname[PETSC_MAX_PATH_LEN];
    321:   PetscMPIInt  rank;

    323:   PetscFunctionBeginUser;
    324:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    325:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
    326:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
    327:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

    329:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank));

    331:   PetscCall(DMGetCoordinateDM(da, &cda));
    332:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    333:   PetscCall(DMDAVecGetArrayRead(cda, coords, &_coords));
    334:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    335:   for (j = sj; j < sj + ny - 1; j++) {
    336:     for (i = si; i < si + nx - 1; i++) {
    337:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
    338:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y)));
    339:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i + 1].x), (double)PetscRealPart(_coords[j + 1][i + 1].y)));
    340:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y)));
    341:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
    342:     }
    343:   }
    344:   PetscCall(DMDAVecRestoreArrayRead(cda, coords, &_coords));

    346:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
    347:   PetscFunctionReturn(PETSC_SUCCESS);
    348: }

    350: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
    351: {
    352:   DM           cda;
    353:   Vec          coords, local_fields;
    354:   DMDACoor2d **_coords;
    355:   FILE        *fp;
    356:   char         fname[PETSC_MAX_PATH_LEN];
    357:   PetscMPIInt  rank;
    358:   PetscInt     si, sj, nx, ny, i, j;
    359:   PetscInt     n_dofs, d;
    360:   PetscScalar *_fields;

    362:   PetscFunctionBeginUser;
    363:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    364:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
    365:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
    366:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

    368:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
    369:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
    370:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
    371:   for (d = 0; d < n_dofs; d++) {
    372:     const char *field_name;
    373:     PetscCall(DMDAGetFieldName(da, d, &field_name));
    374:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
    375:   }
    376:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

    378:   PetscCall(DMGetCoordinateDM(da, &cda));
    379:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    380:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    381:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));

    383:   PetscCall(DMCreateLocalVector(da, &local_fields));
    384:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
    385:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
    386:   PetscCall(VecGetArray(local_fields, &_fields));

    388:   for (j = sj; j < sj + ny; j++) {
    389:     for (i = si; i < si + nx; i++) {
    390:       PetscScalar coord_x, coord_y;
    391:       PetscScalar field_d;

    393:       coord_x = _coords[j][i].x;
    394:       coord_y = _coords[j][i].y;

    396:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
    397:       for (d = 0; d < n_dofs; d++) {
    398:         field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
    399:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d)));
    400:       }
    401:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
    402:     }
    403:   }
    404:   PetscCall(VecRestoreArray(local_fields, &_fields));
    405:   PetscCall(VecDestroy(&local_fields));

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

    409:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
    410:   PetscFunctionReturn(PETSC_SUCCESS);
    411: }

    413: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
    414: {
    415:   DM                       cda;
    416:   Vec                      local_fields;
    417:   FILE                    *fp;
    418:   char                     fname[PETSC_MAX_PATH_LEN];
    419:   PetscMPIInt              rank;
    420:   PetscInt                 si, sj, nx, ny, i, j, p;
    421:   PetscInt                 n_dofs, d;
    422:   GaussPointCoefficients **_coefficients;

    424:   PetscFunctionBeginUser;
    425:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    426:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
    427:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
    428:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

    430:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
    431:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
    432:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
    433:   for (d = 0; d < n_dofs; d++) {
    434:     const char *field_name;
    435:     PetscCall(DMDAGetFieldName(da, d, &field_name));
    436:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
    437:   }
    438:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

    440:   PetscCall(DMGetCoordinateDM(da, &cda));
    441:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));

    443:   PetscCall(DMCreateLocalVector(da, &local_fields));
    444:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
    445:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
    446:   PetscCall(DMDAVecGetArray(da, local_fields, &_coefficients));

    448:   for (j = sj; j < sj + ny; j++) {
    449:     for (i = si; i < si + nx; i++) {
    450:       PetscScalar coord_x, coord_y;

    452:       for (p = 0; p < GAUSS_POINTS; p++) {
    453:         coord_x = _coefficients[j][i].gp_coords[2 * p];
    454:         coord_y = _coefficients[j][i].gp_coords[2 * p + 1];

    456:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));

    458:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e %1.6e", (double)PetscRealPart(_coefficients[j][i].eta[p]), (double)PetscRealPart(_coefficients[j][i].fx[p]), (double)PetscRealPart(_coefficients[j][i].fy[p])));
    459:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
    460:       }
    461:     }
    462:   }
    463:   PetscCall(DMDAVecRestoreArray(da, local_fields, &_coefficients));
    464:   PetscCall(VecDestroy(&local_fields));

    466:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
    467:   PetscFunctionReturn(PETSC_SUCCESS);
    468: }

    470: 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)
    471: {
    472:   PetscInt ij;
    473:   PetscInt r, c, nc;

    475:   nc = u_NPE * u_dof;
    476:   r  = w_dof * wi + wd;
    477:   c  = u_dof * ui + ud;
    478:   ij = r * nc + c;
    479:   return ij;
    480: }

    482: /*
    483:  D = [ 2.eta   0   0   ]
    484:      [   0   2.eta 0   ]
    485:      [   0     0   eta ]

    487:  B = [ d_dx   0   ]
    488:      [  0    d_dy ]
    489:      [ d_dy  d_dx ]
    490: */
    491: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    492: {
    493:   PetscInt    ngp;
    494:   PetscScalar gp_xi[GAUSS_POINTS][2];
    495:   PetscScalar gp_weight[GAUSS_POINTS];
    496:   PetscInt    p, i, j, k;
    497:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    498:   PetscScalar J_p, tildeD[3];
    499:   PetscScalar B[3][U_DOFS * NODES_PER_EL];

    501:   /* define quadrature rule */
    502:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    504:   /* evaluate integral */
    505:   for (p = 0; p < ngp; p++) {
    506:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
    507:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);

    509:     for (i = 0; i < NODES_PER_EL; i++) {
    510:       PetscScalar d_dx_i = GNx_p[0][i];
    511:       PetscScalar d_dy_i = GNx_p[1][i];

    513:       B[0][2 * i]     = d_dx_i;
    514:       B[0][2 * i + 1] = 0.0;
    515:       B[1][2 * i]     = 0.0;
    516:       B[1][2 * i + 1] = d_dy_i;
    517:       B[2][2 * i]     = d_dy_i;
    518:       B[2][2 * i + 1] = d_dx_i;
    519:     }

    521:     tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
    522:     tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
    523:     tildeD[2] = gp_weight[p] * J_p * eta[p];

    525:     /* form Bt tildeD B */
    526:     /*
    527:     Ke_ij = Bt_ik . D_kl . B_lj
    528:           = B_ki . D_kl . B_lj
    529:           = B_ki . D_kk . B_kj
    530:     */
    531:     for (i = 0; i < 8; i++) {
    532:       for (j = 0; j < 8; j++) {
    533:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
    534:           Ke[i + 8 * j] = Ke[i + 8 * j] + B[k][i] * tildeD[k] * B[k][j];
    535:         }
    536:       }
    537:     }
    538:   }
    539: }

    541: static void FormGradientOperatorQ1(PetscScalar Ke[], PetscScalar coords[])
    542: {
    543:   PetscInt    ngp;
    544:   PetscScalar gp_xi[GAUSS_POINTS][2];
    545:   PetscScalar gp_weight[GAUSS_POINTS];
    546:   PetscInt    p, i, j, di;
    547:   PetscScalar Ni_p[NODES_PER_EL];
    548:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    549:   PetscScalar J_p, fac;

    551:   /* define quadrature rule */
    552:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    554:   /* evaluate integral */
    555:   for (p = 0; p < ngp; p++) {
    556:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
    557:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
    558:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
    559:     fac = gp_weight[p] * J_p;

    561:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
    562:       for (di = 0; di < NSD; di++) {     /* u dofs */
    563:         for (j = 0; j < 4; j++) {        /* p nodes, p dofs = 1 (ie no loop) */
    564:           PetscInt IJ;
    565:           IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);

    567:           Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
    568:         }
    569:       }
    570:     }
    571:   }
    572: }

    574: static void FormDivergenceOperatorQ1(PetscScalar De[], PetscScalar coords[])
    575: {
    576:   PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
    577:   PetscInt    i, j;
    578:   PetscInt    nr_g, nc_g;

    580:   PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
    581:   FormGradientOperatorQ1(Ge, coords);

    583:   nr_g = U_DOFS * NODES_PER_EL;
    584:   nc_g = P_DOFS * NODES_PER_EL;

    586:   for (i = 0; i < nr_g; i++) {
    587:     for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
    588:   }
    589: }

    591: static void FormStabilisationOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    592: {
    593:   PetscInt    ngp;
    594:   PetscScalar gp_xi[GAUSS_POINTS][2];
    595:   PetscScalar gp_weight[GAUSS_POINTS];
    596:   PetscInt    p, i, j;
    597:   PetscScalar Ni_p[NODES_PER_EL];
    598:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    599:   PetscScalar J_p, fac, eta_avg;

    601:   /* define quadrature rule */
    602:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    604:   /* evaluate integral */
    605:   for (p = 0; p < ngp; p++) {
    606:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
    607:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
    608:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
    609:     fac = gp_weight[p] * J_p;

    611:     for (i = 0; i < NODES_PER_EL; i++) {
    612:       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] - 0.0625);
    613:     }
    614:   }

    616:   /* scale */
    617:   eta_avg = 0.0;
    618:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    619:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
    620:   fac     = 1.0 / eta_avg;
    621:   for (i = 0; i < NODES_PER_EL; i++) {
    622:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
    623:   }
    624: }

    626: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
    627: {
    628:   PetscInt    ngp;
    629:   PetscScalar gp_xi[GAUSS_POINTS][2];
    630:   PetscScalar gp_weight[GAUSS_POINTS];
    631:   PetscInt    p, i, j;
    632:   PetscScalar Ni_p[NODES_PER_EL];
    633:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    634:   PetscScalar J_p, fac, eta_avg;

    636:   /* define quadrature rule */
    637:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    639:   /* evaluate integral */
    640:   for (p = 0; p < ngp; p++) {
    641:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
    642:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
    643:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
    644:     fac = gp_weight[p] * J_p;

    646:     for (i = 0; i < NODES_PER_EL; i++) {
    647:       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];
    648:     }
    649:   }

    651:   /* scale */
    652:   eta_avg = 0.0;
    653:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    654:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
    655:   fac     = 1.0 / eta_avg;
    656:   for (i = 0; i < NODES_PER_EL; i++) {
    657:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
    658:   }
    659: }

    661: static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
    662: {
    663:   PetscInt    ngp;
    664:   PetscScalar gp_xi[GAUSS_POINTS][2];
    665:   PetscScalar gp_weight[GAUSS_POINTS];
    666:   PetscInt    p, i;
    667:   PetscScalar Ni_p[NODES_PER_EL];
    668:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    669:   PetscScalar J_p, fac;

    671:   /* define quadrature rule */
    672:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    674:   /* evaluate integral */
    675:   for (p = 0; p < ngp; p++) {
    676:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
    677:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
    678:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
    679:     fac = gp_weight[p] * J_p;

    681:     for (i = 0; i < NODES_PER_EL; i++) {
    682:       Fe[NSD * i] += fac * Ni_p[i] * fx[p];
    683:       Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
    684:     }
    685:   }
    686: }

    688: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
    689: {
    690:   PetscFunctionBeginUser;
    691:   /* get coords for the element */
    692:   el_coords[NSD * 0]     = _coords[ej][ei].x;
    693:   el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
    694:   el_coords[NSD * 1]     = _coords[ej + 1][ei].x;
    695:   el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
    696:   el_coords[NSD * 2]     = _coords[ej + 1][ei + 1].x;
    697:   el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
    698:   el_coords[NSD * 3]     = _coords[ej][ei + 1].x;
    699:   el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
    700:   PetscFunctionReturn(PETSC_SUCCESS);
    701: }

    703: static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
    704: {
    705:   DM                       cda;
    706:   Vec                      coords;
    707:   DMDACoor2d             **_coords;
    708:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    709:   MatStencil               p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
    710:   PetscInt                 sex, sey, mx, my;
    711:   PetscInt                 ei, ej;
    712:   PetscScalar              Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
    713:   PetscScalar              Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
    714:   PetscScalar              De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
    715:   PetscScalar              Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
    716:   PetscScalar              el_coords[NODES_PER_EL * NSD];
    717:   Vec                      local_properties;
    718:   GaussPointCoefficients **props;
    719:   PetscScalar             *prop_eta;

    721:   PetscFunctionBeginUser;
    722:   /* setup for coords */
    723:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    724:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    725:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    727:   /* setup for coefficients */
    728:   PetscCall(DMCreateLocalVector(properties_da, &local_properties));
    729:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
    730:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
    731:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

    733:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
    734:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
    735:   for (ej = sey; ej < sey + my; ej++) {
    736:     for (ei = sex; ei < sex + mx; ei++) {
    737:       /* get coords for the element */
    738:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

    740:       /* get coefficients for the element */
    741:       prop_eta = props[ej][ei].eta;

    743:       /* initialise element stiffness matrix */
    744:       PetscCall(PetscMemzero(Ae, sizeof(Ae)));
    745:       PetscCall(PetscMemzero(Ge, sizeof(Ge)));
    746:       PetscCall(PetscMemzero(De, sizeof(De)));
    747:       PetscCall(PetscMemzero(Ce, sizeof(Ce)));

    749:       /* form element stiffness matrix */
    750:       FormStressOperatorQ1(Ae, el_coords, prop_eta);
    751:       FormGradientOperatorQ1(Ge, el_coords);
    752:       FormDivergenceOperatorQ1(De, el_coords);
    753:       FormStabilisationOperatorQ1(Ce, el_coords, prop_eta);

    755:       /* insert element matrix into global matrix */
    756:       PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
    757:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
    758:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
    759:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
    760:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
    761:     }
    762:   }
    763:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    764:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

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

    768:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
    769:   PetscCall(VecDestroy(&local_properties));
    770:   PetscFunctionReturn(PETSC_SUCCESS);
    771: }

    773: static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
    774: {
    775:   DM                       cda;
    776:   Vec                      coords;
    777:   DMDACoor2d             **_coords;
    778:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    779:   MatStencil               p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
    780:   PetscInt                 sex, sey, mx, my;
    781:   PetscInt                 ei, ej;
    782:   PetscScalar              Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
    783:   PetscScalar              Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
    784:   PetscScalar              De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
    785:   PetscScalar              Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
    786:   PetscScalar              el_coords[NODES_PER_EL * NSD];
    787:   Vec                      local_properties;
    788:   GaussPointCoefficients **props;
    789:   PetscScalar             *prop_eta;

    791:   PetscFunctionBeginUser;
    792:   /* setup for coords */
    793:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    794:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    795:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    797:   /* setup for coefficients */
    798:   PetscCall(DMCreateLocalVector(properties_da, &local_properties));
    799:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
    800:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
    801:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

    803:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
    804:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
    805:   for (ej = sey; ej < sey + my; ej++) {
    806:     for (ei = sex; ei < sex + mx; ei++) {
    807:       /* get coords for the element */
    808:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

    810:       /* get coefficients for the element */
    811:       prop_eta = props[ej][ei].eta;

    813:       /* initialise element stiffness matrix */
    814:       PetscCall(PetscMemzero(Ae, sizeof(Ae)));
    815:       PetscCall(PetscMemzero(Ge, sizeof(Ge)));
    816:       PetscCall(PetscMemzero(De, sizeof(De)));
    817:       PetscCall(PetscMemzero(Ce, sizeof(Ce)));

    819:       /* form element stiffness matrix */
    820:       FormStressOperatorQ1(Ae, el_coords, prop_eta);
    821:       FormGradientOperatorQ1(Ge, el_coords);
    822:       FormScaledMassMatrixOperatorQ1(Ce, el_coords, prop_eta);

    824:       /* insert element matrix into global matrix */
    825:       PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
    826:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
    827:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
    828:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
    829:     }
    830:   }
    831:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    832:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

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

    836:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
    837:   PetscCall(VecDestroy(&local_properties));
    838:   PetscFunctionReturn(PETSC_SUCCESS);
    839: }

    841: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
    842: {
    843:   PetscInt n;

    845:   PetscFunctionBeginUser;
    846:   for (n = 0; n < 4; n++) {
    847:     fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].u_dof         = fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].u_dof + Fe_u[2 * n];
    848:     fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].v_dof = fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].v_dof + Fe_u[2 * n + 1];
    849:     fields_F[p_eqn[n].j][p_eqn[n].i].p_dof                 = fields_F[p_eqn[n].j][p_eqn[n].i].p_dof + Fe_p[n];
    850:   }
    851:   PetscFunctionReturn(PETSC_SUCCESS);
    852: }

    854: static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, DM properties_da, Vec properties)
    855: {
    856:   DM                       cda;
    857:   Vec                      coords;
    858:   DMDACoor2d             **_coords;
    859:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    860:   MatStencil               p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
    861:   PetscInt                 sex, sey, mx, my;
    862:   PetscInt                 ei, ej;
    863:   PetscScalar              Fe[NODES_PER_EL * U_DOFS];
    864:   PetscScalar              He[NODES_PER_EL * P_DOFS];
    865:   PetscScalar              el_coords[NODES_PER_EL * NSD];
    866:   Vec                      local_properties;
    867:   GaussPointCoefficients **props;
    868:   PetscScalar             *prop_fx, *prop_fy;
    869:   Vec                      local_F;
    870:   StokesDOF              **ff;

    872:   PetscFunctionBeginUser;
    873:   /* setup for coords */
    874:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    875:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    876:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    878:   /* setup for coefficients */
    879:   PetscCall(DMGetLocalVector(properties_da, &local_properties));
    880:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
    881:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
    882:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

    884:   /* get access to the vector */
    885:   PetscCall(DMGetLocalVector(stokes_da, &local_F));
    886:   PetscCall(VecZeroEntries(local_F));
    887:   PetscCall(DMDAVecGetArray(stokes_da, local_F, &ff));

    889:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
    890:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
    891:   for (ej = sey; ej < sey + my; ej++) {
    892:     for (ei = sex; ei < sex + mx; ei++) {
    893:       /* get coords for the element */
    894:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

    896:       /* get coefficients for the element */
    897:       prop_fx = props[ej][ei].fx;
    898:       prop_fy = props[ej][ei].fy;

    900:       /* initialise element stiffness matrix */
    901:       PetscCall(PetscMemzero(Fe, sizeof(Fe)));
    902:       PetscCall(PetscMemzero(He, sizeof(He)));

    904:       /* form element stiffness matrix */
    905:       FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);

    907:       /* insert element matrix into global matrix */
    908:       PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));

    910:       PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He));
    911:     }
    912:   }

    914:   PetscCall(DMDAVecRestoreArray(stokes_da, local_F, &ff));
    915:   PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
    916:   PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
    917:   PetscCall(DMRestoreLocalVector(stokes_da, &local_F));

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

    921:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
    922:   PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
    923:   PetscFunctionReturn(PETSC_SUCCESS);
    924: }

    926: static PetscErrorCode DMDACreateSolCx(PetscReal eta0, PetscReal eta1, PetscReal xc, PetscInt nz, PetscInt mx, PetscInt my, DM *_da, Vec *_X)
    927: {
    928:   DM           da, cda;
    929:   Vec          X;
    930:   StokesDOF  **_stokes;
    931:   Vec          coords;
    932:   DMDACoor2d **_coords;
    933:   PetscInt     si, sj, ei, ej, i, j;

    935:   PetscFunctionBeginUser;
    936:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, NULL, &da));
    937:   PetscCall(DMSetFromOptions(da));
    938:   PetscCall(DMSetUp(da));
    939:   PetscCall(DMDASetFieldName(da, 0, "anlytic_Vx"));
    940:   PetscCall(DMDASetFieldName(da, 1, "anlytic_Vy"));
    941:   PetscCall(DMDASetFieldName(da, 2, "analytic_P"));

    943:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0., 0.));

    945:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    946:   PetscCall(DMGetCoordinateDM(da, &cda));
    947:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    949:   PetscCall(DMCreateGlobalVector(da, &X));
    950:   PetscCall(DMDAVecGetArray(da, X, &_stokes));

    952:   PetscCall(DMDAGetCorners(da, &si, &sj, 0, &ei, &ej, 0));
    953:   for (j = sj; j < sj + ej; j++) {
    954:     for (i = si; i < si + ei; i++) {
    955:       PetscReal pos[2], pressure, vel[2], total_stress[3], strain_rate[3];

    957:       pos[0] = PetscRealPart(_coords[j][i].x);
    958:       pos[1] = PetscRealPart(_coords[j][i].y);

    960:       evaluate_solCx(pos, eta0, eta1, xc, nz, vel, &pressure, total_stress, strain_rate);

    962:       _stokes[j][i].u_dof = vel[0];
    963:       _stokes[j][i].v_dof = vel[1];
    964:       _stokes[j][i].p_dof = pressure;
    965:     }
    966:   }
    967:   PetscCall(DMDAVecRestoreArray(da, X, &_stokes));
    968:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

    970:   *_da = da;
    971:   *_X  = X;
    972:   PetscFunctionReturn(PETSC_SUCCESS);
    973: }

    975: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields, PetscInt ei, PetscInt ej, StokesDOF nodal_fields[])
    976: {
    977:   PetscFunctionBeginUser;
    978:   /* get the nodal fields */
    979:   nodal_fields[0].u_dof = fields[ej][ei].u_dof;
    980:   nodal_fields[0].v_dof = fields[ej][ei].v_dof;
    981:   nodal_fields[0].p_dof = fields[ej][ei].p_dof;
    982:   nodal_fields[1].u_dof = fields[ej + 1][ei].u_dof;
    983:   nodal_fields[1].v_dof = fields[ej + 1][ei].v_dof;
    984:   nodal_fields[1].p_dof = fields[ej + 1][ei].p_dof;
    985:   nodal_fields[2].u_dof = fields[ej + 1][ei + 1].u_dof;
    986:   nodal_fields[2].v_dof = fields[ej + 1][ei + 1].v_dof;
    987:   nodal_fields[2].p_dof = fields[ej + 1][ei + 1].p_dof;
    988:   nodal_fields[3].u_dof = fields[ej][ei + 1].u_dof;
    989:   nodal_fields[3].v_dof = fields[ej][ei + 1].v_dof;
    990:   nodal_fields[3].p_dof = fields[ej][ei + 1].p_dof;
    991:   PetscFunctionReturn(PETSC_SUCCESS);
    992: }

    994: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da, Vec X, Vec X_analytic)
    995: {
    996:   DM           cda;
    997:   Vec          coords, X_analytic_local, X_local;
    998:   DMDACoor2d **_coords;
    999:   PetscInt     sex, sey, mx, my;
    1000:   PetscInt     ei, ej;
    1001:   PetscScalar  el_coords[NODES_PER_EL * NSD];
    1002:   StokesDOF  **stokes_analytic, **stokes;
    1003:   StokesDOF    stokes_analytic_e[4], stokes_e[4];

    1005:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    1006:   PetscScalar Ni_p[NODES_PER_EL];
    1007:   PetscInt    ngp;
    1008:   PetscScalar gp_xi[GAUSS_POINTS][2];
    1009:   PetscScalar gp_weight[GAUSS_POINTS];
    1010:   PetscInt    p, i;
    1011:   PetscScalar J_p, fac;
    1012:   PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
    1013:   PetscInt    M;
    1014:   PetscReal   xymin[2], xymax[2];

    1016:   PetscFunctionBeginUser;
    1017:   /* define quadrature rule */
    1018:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    1020:   /* setup for coords */
    1021:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
    1022:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
    1023:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    1025:   /* setup for analytic */
    1026:   PetscCall(DMCreateLocalVector(stokes_da, &X_analytic_local));
    1027:   PetscCall(DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
    1028:   PetscCall(DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
    1029:   PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));

    1031:   /* setup for solution */
    1032:   PetscCall(DMCreateLocalVector(stokes_da, &X_local));
    1033:   PetscCall(DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local));
    1034:   PetscCall(DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local));
    1035:   PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));

    1037:   PetscCall(DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    1038:   PetscCall(DMGetBoundingBox(stokes_da, xymin, xymax));

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

    1042:   tp_L2 = tu_L2 = tu_H1 = 0.0;

    1044:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
    1045:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
    1046:   for (ej = sey; ej < sey + my; ej++) {
    1047:     for (ei = sex; ei < sex + mx; ei++) {
    1048:       /* get coords for the element */
    1049:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
    1050:       PetscCall(StokesDAGetNodalFields(stokes, ei, ej, stokes_e));
    1051:       PetscCall(StokesDAGetNodalFields(stokes_analytic, ei, ej, stokes_analytic_e));

    1053:       /* evaluate integral */
    1054:       p_e_L2 = 0.0;
    1055:       u_e_L2 = 0.0;
    1056:       u_e_H1 = 0.0;
    1057:       for (p = 0; p < ngp; p++) {
    1058:         ConstructQ12D_Ni(gp_xi[p], Ni_p);
    1059:         ConstructQ12D_GNi(gp_xi[p], GNi_p);
    1060:         ConstructQ12D_GNx(GNi_p, GNx_p, el_coords, &J_p);
    1061:         fac = gp_weight[p] * J_p;

    1063:         for (i = 0; i < NODES_PER_EL; i++) {
    1064:           PetscScalar u_error, v_error;

    1066:           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);

    1068:           u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
    1069:           v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
    1070:           u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error);

    1072:           u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error     /* du/dx */
    1073:                                    + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error   /* du/dy */
    1074:                                    + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error   /* dv/dx */
    1075:                                    + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error); /* dv/dy */
    1076:         }
    1077:       }

    1079:       tp_L2 += p_e_L2;
    1080:       tu_L2 += u_e_L2;
    1081:       tu_H1 += u_e_H1;
    1082:     }
    1083:   }
    1084:   PetscCallMPI(MPIU_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    1085:   PetscCallMPI(MPIU_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    1086:   PetscCallMPI(MPIU_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
    1087:   p_L2 = PetscSqrtScalar(p_L2);
    1088:   u_L2 = PetscSqrtScalar(u_L2);
    1089:   u_H1 = PetscSqrtScalar(u_H1);

    1091:   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)));

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

    1095:   PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
    1096:   PetscCall(VecDestroy(&X_analytic_local));
    1097:   PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
    1098:   PetscCall(VecDestroy(&X_local));
    1099:   PetscFunctionReturn(PETSC_SUCCESS);
    1100: }

    1102: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx, PetscInt my)
    1103: {
    1104:   DM                       da_Stokes, da_prop;
    1105:   PetscInt                 u_dof, p_dof, dof, stencil_width;
    1106:   Mat                      A, B;
    1107:   DM                       prop_cda, vel_cda;
    1108:   Vec                      prop_coords, vel_coords;
    1109:   PetscInt                 si, sj, nx, ny, i, j, p;
    1110:   Vec                      f, X;
    1111:   PetscInt                 prop_dof, prop_stencil_width;
    1112:   Vec                      properties, l_properties;
    1113:   PetscReal                dx, dy;
    1114:   PetscInt                 M, N;
    1115:   DMDACoor2d             **_prop_coords, **_vel_coords;
    1116:   GaussPointCoefficients **element_props;
    1117:   PetscInt                 its;
    1118:   KSP                      ksp_S;
    1119:   PetscInt                 coefficient_structure = 0;
    1120:   PetscInt                 cpu_x, cpu_y, *lx = NULL, *ly = NULL;
    1121:   PetscBool                use_gp_coords = PETSC_FALSE, set, output_gnuplot = PETSC_FALSE, glvis = PETSC_FALSE, change = PETSC_FALSE;
    1122:   char                     filename[PETSC_MAX_PATH_LEN];

    1124:   PetscFunctionBeginUser;
    1125:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-gnuplot", &output_gnuplot, NULL));
    1126:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-glvis", &glvis, NULL));
    1127:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-change", &change, NULL));

    1129:   /* Generate the da for velocity and pressure */
    1130:   /*
    1131:   We use Q1 elements for the temperature.
    1132:   FEM has a 9-point stencil (BOX) or connectivity pattern
    1133:   Num nodes in each direction is mx+1, my+1
    1134:   */
    1135:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
    1136:   p_dof         = P_DOFS; /* p - pressure */
    1137:   dof           = u_dof + p_dof;
    1138:   stencil_width = 1;
    1139:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &da_Stokes));

    1141:   PetscCall(DMSetMatType(da_Stokes, MATAIJ));
    1142:   PetscCall(DMSetFromOptions(da_Stokes));
    1143:   PetscCall(DMSetUp(da_Stokes));
    1144:   PetscCall(DMDASetFieldName(da_Stokes, 0, "Vx"));
    1145:   PetscCall(DMDASetFieldName(da_Stokes, 1, "Vy"));
    1146:   PetscCall(DMDASetFieldName(da_Stokes, 2, "P"));

    1148:   /* unit box [0,1] x [0,1] */
    1149:   PetscCall(DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0., 0.));

    1151:   /* Generate element properties, we will assume all material properties are constant over the element */
    1152:   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!!  */
    1153:   PetscCall(DMDAGetInfo(da_Stokes, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
    1154:   PetscCall(DMDAGetElementOwnershipRanges2d(da_Stokes, &lx, &ly));

    1156:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
    1157:   prop_stencil_width = 0;
    1158:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, cpu_x, cpu_y, prop_dof, prop_stencil_width, lx, ly, &da_prop));
    1159:   PetscCall(DMSetFromOptions(da_prop));
    1160:   PetscCall(DMSetUp(da_prop));
    1161:   PetscCall(PetscFree(lx));
    1162:   PetscCall(PetscFree(ly));

    1164:   /* define centroid positions */
    1165:   PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    1166:   dx = 1.0 / (PetscReal)M;
    1167:   dy = 1.0 / (PetscReal)N;

    1169:   PetscCall(DMDASetUniformCoordinates(da_prop, 0.0 + 0.5 * dx, 1.0 - 0.5 * dx, 0.0 + 0.5 * dy, 1.0 - 0.5 * dy, 0., 0));

    1171:   /* define coefficients */
    1172:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));

    1174:   PetscCall(DMCreateGlobalVector(da_prop, &properties));
    1175:   PetscCall(DMCreateLocalVector(da_prop, &l_properties));
    1176:   PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));

    1178:   PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
    1179:   PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
    1180:   PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));

    1182:   PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));

    1184:   PetscCall(DMGetCoordinateDM(da_Stokes, &vel_cda));
    1185:   PetscCall(DMGetCoordinatesLocal(da_Stokes, &vel_coords));
    1186:   PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));

    1188:   /* interpolate the coordinates */
    1189:   for (j = sj; j < sj + ny; j++) {
    1190:     for (i = si; i < si + nx; i++) {
    1191:       PetscInt    ngp;
    1192:       PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
    1193:       PetscScalar el_coords[8];

    1195:       PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));
    1196:       ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    1198:       for (p = 0; p < GAUSS_POINTS; p++) {
    1199:         PetscScalar gp_x, gp_y;
    1200:         PetscInt    n;
    1201:         PetscScalar xi_p[2], Ni_p[4];

    1203:         xi_p[0] = gp_xi[p][0];
    1204:         xi_p[1] = gp_xi[p][1];
    1205:         ConstructQ12D_Ni(xi_p, Ni_p);

    1207:         gp_x = 0.0;
    1208:         gp_y = 0.0;
    1209:         for (n = 0; n < NODES_PER_EL; n++) {
    1210:           gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
    1211:           gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
    1212:         }
    1213:         element_props[j][i].gp_coords[2 * p]     = gp_x;
    1214:         element_props[j][i].gp_coords[2 * p + 1] = gp_y;
    1215:       }
    1216:     }
    1217:   }

    1219:   /* define the coefficients */
    1220:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, NULL));

    1222:   for (j = sj; j < sj + ny; j++) {
    1223:     for (i = si; i < si + nx; i++) {
    1224:       PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
    1225:       PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
    1226:       PetscReal coord_x, coord_y;

    1228:       if (coefficient_structure == 0) {
    1229:         PetscReal opts_eta0, opts_eta1, opts_xc;
    1230:         PetscInt  opts_nz;

    1232:         opts_eta0 = 1.0;
    1233:         opts_eta1 = 1.0;
    1234:         opts_xc   = 0.5;
    1235:         opts_nz   = 1;

    1237:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
    1238:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
    1239:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
    1240:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));

    1242:         for (p = 0; p < GAUSS_POINTS; p++) {
    1243:           coord_x = centroid_x;
    1244:           coord_y = centroid_y;
    1245:           if (use_gp_coords) {
    1246:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
    1247:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
    1248:           }

    1250:           element_props[j][i].eta[p] = opts_eta0;
    1251:           if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;

    1253:           element_props[j][i].fx[p] = 0.0;
    1254:           element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
    1255:         }
    1256:       } else if (coefficient_structure == 1) { /* square sinker */
    1257:         PetscReal opts_eta0, opts_eta1, opts_dx, opts_dy;

    1259:         opts_eta0 = 1.0;
    1260:         opts_eta1 = 1.0;
    1261:         opts_dx   = 0.50;
    1262:         opts_dy   = 0.50;

    1264:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
    1265:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
    1266:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
    1267:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));

    1269:         for (p = 0; p < GAUSS_POINTS; p++) {
    1270:           coord_x = centroid_x;
    1271:           coord_y = centroid_y;
    1272:           if (use_gp_coords) {
    1273:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
    1274:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
    1275:           }

    1277:           element_props[j][i].eta[p] = opts_eta0;
    1278:           element_props[j][i].fx[p]  = 0.0;
    1279:           element_props[j][i].fy[p]  = 0.0;

    1281:           if ((coord_x > -0.5 * opts_dx + 0.5) && (coord_x < 0.5 * opts_dx + 0.5)) {
    1282:             if ((coord_y > -0.5 * opts_dy + 0.5) && (coord_y < 0.5 * opts_dy + 0.5)) {
    1283:               element_props[j][i].eta[p] = opts_eta1;
    1284:               element_props[j][i].fx[p]  = 0.0;
    1285:               element_props[j][i].fy[p]  = -1.0;
    1286:             }
    1287:           }
    1288:         }
    1289:       } else if (coefficient_structure == 2) { /* circular sinker */
    1290:         PetscReal opts_eta0, opts_eta1, opts_r, radius2;

    1292:         opts_eta0 = 1.0;
    1293:         opts_eta1 = 1.0;
    1294:         opts_r    = 0.25;

    1296:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
    1297:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
    1298:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));

    1300:         for (p = 0; p < GAUSS_POINTS; p++) {
    1301:           coord_x = centroid_x;
    1302:           coord_y = centroid_y;
    1303:           if (use_gp_coords) {
    1304:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
    1305:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
    1306:           }

    1308:           element_props[j][i].eta[p] = opts_eta0;
    1309:           element_props[j][i].fx[p]  = 0.0;
    1310:           element_props[j][i].fy[p]  = 0.0;

    1312:           radius2 = (coord_x - 0.5) * (coord_x - 0.5) + (coord_y - 0.5) * (coord_y - 0.5);
    1313:           if (radius2 < opts_r * opts_r) {
    1314:             element_props[j][i].eta[p] = opts_eta1;
    1315:             element_props[j][i].fx[p]  = 0.0;
    1316:             element_props[j][i].fy[p]  = -1.0;
    1317:           }
    1318:         }
    1319:       } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
    1320:         PetscReal opts_eta0, opts_eta1, opts_r, opts_dx, opts_dy, opts_c0x, opts_c0y, opts_s0x, opts_s0y, opts_phi, radius2;

    1322:         opts_eta0 = 1.0;
    1323:         opts_eta1 = 1.0;
    1324:         opts_r    = 0.25;
    1325:         opts_c0x  = 0.35; /* circle center */
    1326:         opts_c0y  = 0.35;
    1327:         opts_s0x  = 0.7; /* square center */
    1328:         opts_s0y  = 0.7;
    1329:         opts_dx   = 0.25;
    1330:         opts_dy   = 0.25;
    1331:         opts_phi  = 25;

    1333:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
    1334:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
    1335:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
    1336:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0x", &opts_c0x, NULL));
    1337:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0y", &opts_c0y, NULL));
    1338:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0x", &opts_s0x, NULL));
    1339:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0y", &opts_s0y, NULL));
    1340:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
    1341:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
    1342:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_phi", &opts_phi, NULL));
    1343:         opts_phi *= PETSC_PI / 180;

    1345:         for (p = 0; p < GAUSS_POINTS; p++) {
    1346:           coord_x = centroid_x;
    1347:           coord_y = centroid_y;
    1348:           if (use_gp_coords) {
    1349:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
    1350:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
    1351:           }

    1353:           element_props[j][i].eta[p] = opts_eta0;
    1354:           element_props[j][i].fx[p]  = 0.0;
    1355:           element_props[j][i].fy[p]  = -0.2;

    1357:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
    1358:           if (radius2 < opts_r * opts_r || (PetscAbs(+(coord_x - opts_s0x) * PetscCosReal(opts_phi) + (coord_y - opts_s0y) * PetscSinReal(opts_phi)) < opts_dx / 2 && PetscAbs(-(coord_x - opts_s0x) * PetscSinReal(opts_phi) + (coord_y - opts_s0y) * PetscCosReal(opts_phi)) < opts_dy / 2)) {
    1359:             element_props[j][i].eta[p] = opts_eta1;
    1360:             element_props[j][i].fx[p]  = 0.0;
    1361:             element_props[j][i].fy[p]  = -1.0;
    1362:           }
    1363:         }
    1364:       } else if (coefficient_structure == 4) { /* subdomain jump */
    1365:         PetscReal opts_mag, opts_eta0;
    1366:         PetscInt  opts_nz, px, py;
    1367:         PetscBool jump;

    1369:         opts_mag  = 1.0;
    1370:         opts_eta0 = 1.0;
    1371:         opts_nz   = 1;

    1373:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL));
    1374:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL));
    1375:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-jump_nz", &opts_nz, NULL));
    1376:         PetscCall(DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    1377:         if (px % 2) {
    1378:           jump = (PetscBool)(PetscGlobalRank % 2);
    1379:         } else {
    1380:           jump = (PetscBool)((PetscGlobalRank / px) % 2 ? PetscGlobalRank % 2 : !(PetscGlobalRank % 2));
    1381:         }
    1382:         for (p = 0; p < GAUSS_POINTS; p++) {
    1383:           coord_x = centroid_x;
    1384:           coord_y = centroid_y;
    1385:           if (use_gp_coords) {
    1386:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
    1387:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
    1388:           }

    1390:           element_props[j][i].eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
    1391:           element_props[j][i].fx[p]  = 0.0;
    1392:           element_props[j][i].fy[p]  = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
    1393:         }
    1394:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
    1395:     }
    1396:   }
    1397:   PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));

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

    1401:   PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
    1402:   PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
    1403:   PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));

    1405:   if (output_gnuplot) {
    1406:     PetscCall(DMDACoordViewGnuplot2d(da_Stokes, "mesh"));
    1407:     PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for Stokes eqn.", "properties"));
    1408:   }

    1410:   if (glvis) {
    1411:     Vec         glv_prop, etaf;
    1412:     PetscViewer view;
    1413:     PetscInt    dim;
    1414:     const char *fec = {"FiniteElementCollection: L2_2D_P1"};

    1416:     PetscCall(DMGetDimension(da_Stokes, &dim));
    1417:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, GAUSS_POINTS * mx * mx, &etaf));
    1418:     PetscCall(PetscObjectSetName((PetscObject)etaf, "viscosity"));
    1419:     PetscCall(PetscViewerGLVisOpen(PETSC_COMM_WORLD, PETSC_VIEWER_GLVIS_SOCKET, NULL, PETSC_DECIDE, &view));
    1420:     PetscCall(PetscViewerGLVisSetFields(view, 1, &fec, &dim, glvis_extract_eta, (PetscObject *)&etaf, da_prop, NULL));
    1421:     PetscCall(DMGetLocalVector(da_prop, &glv_prop));
    1422:     PetscCall(DMGlobalToLocalBegin(da_prop, properties, INSERT_VALUES, glv_prop));
    1423:     PetscCall(DMGlobalToLocalEnd(da_prop, properties, INSERT_VALUES, glv_prop));
    1424:     PetscCall(VecSetDM(etaf, da_Stokes));
    1425:     PetscCall(VecView(glv_prop, view));
    1426:     PetscCall(PetscViewerDestroy(&view));
    1427:     PetscCall(DMRestoreLocalVector(da_prop, &glv_prop));
    1428:     PetscCall(VecDestroy(&etaf));
    1429:   }

    1431:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
    1432:   PetscCall(DMCreateMatrix(da_Stokes, &A));
    1433:   PetscCall(DMCreateMatrix(da_Stokes, &B));
    1434:   PetscCall(DMCreateGlobalVector(da_Stokes, &f));
    1435:   PetscCall(DMCreateGlobalVector(da_Stokes, &X));

    1437:   /* assemble A11 */
    1438:   PetscCall(MatZeroEntries(A));
    1439:   PetscCall(MatZeroEntries(B));
    1440:   PetscCall(VecZeroEntries(f));

    1442:   PetscCall(AssembleA_Stokes(A, da_Stokes, da_prop, properties));
    1443:   PetscCall(AssembleA_PCStokes(B, da_Stokes, da_prop, properties));
    1444:   /* build force vector */
    1445:   PetscCall(AssembleF_Stokes(f, da_Stokes, da_prop, properties));

    1447:   PetscCall(DMDABCApplyFreeSlip(da_Stokes, A, f));
    1448:   PetscCall(DMDABCApplyFreeSlip(da_Stokes, B, NULL));

    1450:   /* SOLVE */
    1451:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_S));
    1452:   PetscCall(KSPSetOptionsPrefix(ksp_S, "stokes_"));
    1453:   PetscCall(KSPSetDM(ksp_S, da_Stokes));
    1454:   PetscCall(KSPSetDMActive(ksp_S, PETSC_FALSE));
    1455:   PetscCall(KSPSetOperators(ksp_S, A, B));
    1456:   PetscCall(KSPSetFromOptions(ksp_S));
    1457:   {
    1458:     PC             pc;
    1459:     const PetscInt ufields[] = {0, 1}, pfields[1] = {2};

    1461:     PetscCall(KSPGetPC(ksp_S, &pc));
    1462:     PetscCall(PCFieldSplitSetBlockSize(pc, 3));
    1463:     PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
    1464:     PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
    1465:   }

    1467:   {
    1468:     PC        pc;
    1469:     PetscBool same = PETSC_FALSE;
    1470:     PetscCall(KSPGetPC(ksp_S, &pc));
    1471:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same));
    1472:     if (same) {
    1473:       PetscBool usedivmat = PETSC_FALSE;
    1474:       PetscCall(KSPSetOperators(ksp_S, A, A));

    1476:       PetscCall(PetscOptionsGetBool(NULL, NULL, "-stokes_pc_bddc_use_divergence", &usedivmat, NULL));
    1477:       if (usedivmat) {
    1478:         IS      *fields, vel;
    1479:         PetscInt i, nf;

    1481:         PetscCall(DMCreateFieldDecomposition(da_Stokes, &nf, NULL, &fields, NULL));
    1482:         PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, fields, &vel));
    1483:         PetscCall(MatZeroRowsIS(B, fields[2], 1.0, NULL, NULL)); /* we put 1.0 on the diagonal to pick the pressure average too */
    1484:         PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B));
    1485:         PetscCall(MatZeroRowsIS(B, vel, 0.0, NULL, NULL));
    1486:         PetscCall(ISDestroy(&vel));
    1487:         PetscCall(PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL));
    1488:         for (i = 0; i < nf; i++) PetscCall(ISDestroy(&fields[i]));
    1489:         PetscCall(PetscFree(fields));
    1490:       }
    1491:     }
    1492:   }

    1494:   {
    1495:     PC        pc_S;
    1496:     KSP      *sub_ksp, ksp_U;
    1497:     PetscInt  nsplits;
    1498:     DM        da_U;
    1499:     PetscBool is_pcfs;

    1501:     PetscCall(KSPSetUp(ksp_S));
    1502:     PetscCall(KSPGetPC(ksp_S, &pc_S));

    1504:     is_pcfs = PETSC_FALSE;
    1505:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_S, PCFIELDSPLIT, &is_pcfs));

    1507:     if (is_pcfs) {
    1508:       PetscCall(PCFieldSplitGetSubKSP(pc_S, &nsplits, &sub_ksp));
    1509:       ksp_U = sub_ksp[0];
    1510:       PetscCall(PetscFree(sub_ksp));

    1512:       if (nsplits == 2) {
    1513:         PetscCall(DMDACreateCompatibleDMDA(da_Stokes, 2, &da_U));

    1515:         PetscCall(KSPSetDM(ksp_U, da_U));
    1516:         PetscCall(KSPSetDMActive(ksp_U, PETSC_FALSE));
    1517:         PetscCall(DMDestroy(&da_U));
    1518:         if (change) {
    1519:           const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
    1520:           PetscCall(PetscOptionsInsertString(NULL, opt));
    1521:           PetscCall(KSPSetFromOptions(ksp_U));
    1522:         }
    1523:         PetscCall(KSPSetUp(ksp_U));
    1524:       }
    1525:     }
    1526:   }

    1528:   PetscCall(KSPSolve(ksp_S, f, X));

    1530:   PetscCall(PetscOptionsGetString(NULL, NULL, "-o", filename, sizeof(filename), &set));
    1531:   if (set) {
    1532:     char       *ext;
    1533:     PetscViewer viewer;
    1534:     PetscBool   flg;
    1535:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
    1536:     PetscCall(PetscStrrchr(filename, '.', &ext));
    1537:     PetscCall(PetscStrcmp("vts", ext, &flg));
    1538:     if (flg) {
    1539:       PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
    1540:     } else {
    1541:       PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
    1542:     }
    1543:     PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
    1544:     PetscCall(PetscViewerFileSetName(viewer, filename));
    1545:     PetscCall(VecView(X, viewer));
    1546:     PetscCall(PetscViewerDestroy(&viewer));
    1547:   }
    1548:   if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes, X, "Velocity solution for Stokes eqn.", "X"));

    1550:   if (glvis) {
    1551:     PetscViewer view;

    1553:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &view));
    1554:     PetscCall(PetscViewerSetType(view, PETSCVIEWERGLVIS));
    1555:     PetscCall(VecView(X, view));
    1556:     PetscCall(PetscViewerDestroy(&view));
    1557:   }

    1559:   PetscCall(KSPGetIterationNumber(ksp_S, &its));

    1561:   if (coefficient_structure == 0) {
    1562:     PetscReal opts_eta0, opts_eta1, opts_xc;
    1563:     PetscInt  opts_nz, N;
    1564:     DM        da_Stokes_analytic;
    1565:     Vec       X_analytic;
    1566:     PetscReal nrm1[3], nrm2[3], nrmI[3];

    1568:     opts_eta0 = 1.0;
    1569:     opts_eta1 = 1.0;
    1570:     opts_xc   = 0.5;
    1571:     opts_nz   = 1;

    1573:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
    1574:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
    1575:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
    1576:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));

    1578:     PetscCall(DMDACreateSolCx(opts_eta0, opts_eta1, opts_xc, opts_nz, mx, my, &da_Stokes_analytic, &X_analytic));
    1579:     if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes_analytic, X_analytic, "Analytic solution for Stokes eqn.", "X_analytic"));
    1580:     PetscCall(DMDAIntegrateErrors(da_Stokes_analytic, X, X_analytic));

    1582:     PetscCall(VecAXPY(X_analytic, -1.0, X));
    1583:     PetscCall(VecGetSize(X_analytic, &N));
    1584:     N = N / 3;

    1586:     PetscCall(VecStrideNorm(X_analytic, 0, NORM_1, &nrm1[0]));
    1587:     PetscCall(VecStrideNorm(X_analytic, 0, NORM_2, &nrm2[0]));
    1588:     PetscCall(VecStrideNorm(X_analytic, 0, NORM_INFINITY, &nrmI[0]));

    1590:     PetscCall(VecStrideNorm(X_analytic, 1, NORM_1, &nrm1[1]));
    1591:     PetscCall(VecStrideNorm(X_analytic, 1, NORM_2, &nrm2[1]));
    1592:     PetscCall(VecStrideNorm(X_analytic, 1, NORM_INFINITY, &nrmI[1]));

    1594:     PetscCall(VecStrideNorm(X_analytic, 2, NORM_1, &nrm1[2]));
    1595:     PetscCall(VecStrideNorm(X_analytic, 2, NORM_2, &nrm2[2]));
    1596:     PetscCall(VecStrideNorm(X_analytic, 2, NORM_INFINITY, &nrmI[2]));

    1598:     PetscCall(DMDestroy(&da_Stokes_analytic));
    1599:     PetscCall(VecDestroy(&X_analytic));
    1600:   }

    1602:   PetscCall(KSPDestroy(&ksp_S));
    1603:   PetscCall(VecDestroy(&X));
    1604:   PetscCall(VecDestroy(&f));
    1605:   PetscCall(MatDestroy(&A));
    1606:   PetscCall(MatDestroy(&B));

    1608:   PetscCall(DMDestroy(&da_Stokes));
    1609:   PetscCall(DMDestroy(&da_prop));

    1611:   PetscCall(VecDestroy(&properties));
    1612:   PetscCall(VecDestroy(&l_properties));
    1613:   PetscFunctionReturn(PETSC_SUCCESS);
    1614: }

    1616: int main(int argc, char **args)
    1617: {
    1618:   PetscInt mx, my;

    1620:   PetscFunctionBeginUser;
    1621:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
    1622:   mx = my = 10;
    1623:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
    1624:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
    1625:   PetscCall(solve_stokes_2d_coupled(mx, my));
    1626:   PetscCall(PetscFinalize());
    1627:   return 0;
    1628: }

    1630: /* -------------------------- helpers for boundary conditions -------------------------------- */
    1631: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
    1632: {
    1633:   DM                     cda;
    1634:   Vec                    coords;
    1635:   PetscInt               si, sj, nx, ny, i, j;
    1636:   PetscInt               M, N;
    1637:   DMDACoor2d           **_coords;
    1638:   const PetscInt        *g_idx;
    1639:   PetscInt              *bc_global_ids;
    1640:   PetscScalar           *bc_vals;
    1641:   PetscInt               nbcs;
    1642:   PetscInt               n_dofs;
    1643:   ISLocalToGlobalMapping ltogm;

    1645:   PetscFunctionBeginUser;
    1646:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1647:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1649:   PetscCall(DMGetCoordinateDM(da, &cda));
    1650:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1651:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1652:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1653:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1655:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
    1656:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

    1658:   /* init the entries to -1 so VecSetValues will ignore them */
    1659:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

    1661:   i = nx - 1;
    1662:   for (j = 0; j < ny; j++) {
    1663:     PetscInt local_id;

    1665:     local_id = i + j * nx;

    1667:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

    1669:     bc_vals[j] = 0.0;
    1670:   }
    1671:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1672:   nbcs = 0;
    1673:   if ((si + nx) == (M)) nbcs = ny;

    1675:   if (b) {
    1676:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1677:     PetscCall(VecAssemblyBegin(b));
    1678:     PetscCall(VecAssemblyEnd(b));
    1679:   }
    1680:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1682:   PetscCall(PetscFree(bc_vals));
    1683:   PetscCall(PetscFree(bc_global_ids));

    1685:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1686:   PetscFunctionReturn(PETSC_SUCCESS);
    1687: }

    1689: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
    1690: {
    1691:   DM                     cda;
    1692:   Vec                    coords;
    1693:   PetscInt               si, sj, nx, ny, i, j;
    1694:   PetscInt               M, N;
    1695:   DMDACoor2d           **_coords;
    1696:   const PetscInt        *g_idx;
    1697:   PetscInt              *bc_global_ids;
    1698:   PetscScalar           *bc_vals;
    1699:   PetscInt               nbcs;
    1700:   PetscInt               n_dofs;
    1701:   ISLocalToGlobalMapping ltogm;

    1703:   PetscFunctionBeginUser;
    1704:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1705:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1707:   PetscCall(DMGetCoordinateDM(da, &cda));
    1708:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1709:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1710:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1711:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1713:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
    1714:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

    1716:   /* init the entries to -1 so VecSetValues will ignore them */
    1717:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

    1719:   i = 0;
    1720:   for (j = 0; j < ny; j++) {
    1721:     PetscInt local_id;

    1723:     local_id = i + j * nx;

    1725:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

    1727:     bc_vals[j] = 0.0;
    1728:   }
    1729:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1730:   nbcs = 0;
    1731:   if (si == 0) nbcs = ny;

    1733:   if (b) {
    1734:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1735:     PetscCall(VecAssemblyBegin(b));
    1736:     PetscCall(VecAssemblyEnd(b));
    1737:   }

    1739:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1741:   PetscCall(PetscFree(bc_vals));
    1742:   PetscCall(PetscFree(bc_global_ids));

    1744:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1745:   PetscFunctionReturn(PETSC_SUCCESS);
    1746: }

    1748: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
    1749: {
    1750:   DM                     cda;
    1751:   Vec                    coords;
    1752:   PetscInt               si, sj, nx, ny, i, j;
    1753:   PetscInt               M, N;
    1754:   DMDACoor2d           **_coords;
    1755:   const PetscInt        *g_idx;
    1756:   PetscInt              *bc_global_ids;
    1757:   PetscScalar           *bc_vals;
    1758:   PetscInt               nbcs;
    1759:   PetscInt               n_dofs;
    1760:   ISLocalToGlobalMapping ltogm;

    1762:   PetscFunctionBeginUser;
    1763:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1764:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1766:   PetscCall(DMGetCoordinateDM(da, &cda));
    1767:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1768:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1769:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1770:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1772:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
    1773:   PetscCall(PetscMalloc1(nx, &bc_vals));

    1775:   /* init the entries to -1 so VecSetValues will ignore them */
    1776:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

    1778:   j = ny - 1;
    1779:   for (i = 0; i < nx; i++) {
    1780:     PetscInt local_id;

    1782:     local_id = i + j * nx;

    1784:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

    1786:     bc_vals[i] = 0.0;
    1787:   }
    1788:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1789:   nbcs = 0;
    1790:   if ((sj + ny) == (N)) nbcs = nx;

    1792:   if (b) {
    1793:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1794:     PetscCall(VecAssemblyBegin(b));
    1795:     PetscCall(VecAssemblyEnd(b));
    1796:   }
    1797:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1799:   PetscCall(PetscFree(bc_vals));
    1800:   PetscCall(PetscFree(bc_global_ids));

    1802:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1803:   PetscFunctionReturn(PETSC_SUCCESS);
    1804: }

    1806: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
    1807: {
    1808:   DM                     cda;
    1809:   Vec                    coords;
    1810:   PetscInt               si, sj, nx, ny, i, j;
    1811:   PetscInt               M, N;
    1812:   DMDACoor2d           **_coords;
    1813:   const PetscInt        *g_idx;
    1814:   PetscInt              *bc_global_ids;
    1815:   PetscScalar           *bc_vals;
    1816:   PetscInt               nbcs;
    1817:   PetscInt               n_dofs;
    1818:   ISLocalToGlobalMapping ltogm;

    1820:   PetscFunctionBeginUser;
    1821:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1822:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1824:   PetscCall(DMGetCoordinateDM(da, &cda));
    1825:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1826:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1827:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1828:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1830:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
    1831:   PetscCall(PetscMalloc1(nx, &bc_vals));

    1833:   /* init the entries to -1 so VecSetValues will ignore them */
    1834:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

    1836:   j = 0;
    1837:   for (i = 0; i < nx; i++) {
    1838:     PetscInt local_id;

    1840:     local_id = i + j * nx;

    1842:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

    1844:     bc_vals[i] = 0.0;
    1845:   }
    1846:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1847:   nbcs = 0;
    1848:   if (sj == 0) nbcs = nx;

    1850:   if (b) {
    1851:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1852:     PetscCall(VecAssemblyBegin(b));
    1853:     PetscCall(VecAssemblyEnd(b));
    1854:   }
    1855:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1857:   PetscCall(PetscFree(bc_vals));
    1858:   PetscCall(PetscFree(bc_global_ids));

    1860:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1861:   PetscFunctionReturn(PETSC_SUCCESS);
    1862: }

    1864: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
    1865: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes, Mat A, Vec f)
    1866: {
    1867:   PetscFunctionBeginUser;
    1868:   PetscCall(BCApplyZero_NORTH(da_Stokes, 1, A, f));
    1869:   PetscCall(BCApplyZero_EAST(da_Stokes, 0, A, f));
    1870:   PetscCall(BCApplyZero_SOUTH(da_Stokes, 1, A, f));
    1871:   PetscCall(BCApplyZero_WEST(da_Stokes, 0, A, f));
    1872:   PetscFunctionReturn(PETSC_SUCCESS);
    1873: }

    1875: /*TEST

    1877:    build:
    1878:       requires: !complex !single

    1880:    test:
    1881:       args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short

    1883:    testset:
    1884:       args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_fieldsplit_u_ksp_type preonly -stokes_fieldsplit_u_pc_type lu -stokes_fieldsplit_p_ksp_type preonly -stokes_fieldsplit_p_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short
    1885:       test:
    1886:          suffix: 2
    1887:          args:
    1888:          output_file: output/ex43_1.out
    1889:       test:
    1890:          requires: mumps
    1891:          suffix: 2_mumps
    1892:          args: -change true -stokes_ksp_view
    1893:          output_file: output/ex43_2_mumps.out
    1894:          # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
    1895:          filter:  grep -E -v "(INFOG\([^7]|INFO\(|\[0\])" | grep -v "            total: nonzeros="

    1897:    test:
    1898:       suffix: 3
    1899:       nsize: 4
    1900:       args: -stokes_pc_use_amat -stokes_ksp_type gcr -stokes_ksp_gcr_restart 60 -stokes_ksp_norm_type unpreconditioned -stokes_ksp_rtol 1.e-2 -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view

    1902:    test:
    1903:       suffix: 4
    1904:       nsize: 4
    1905:       args: -stokes_ksp_type pipegcr -stokes_ksp_pipegcr_mmax 60 -stokes_ksp_pipegcr_unroll_w 1 -stokes_ksp_norm_type natural -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view

    1907:    test:
    1908:       suffix: 5
    1909:       nsize: 4
    1910:       args: -stokes_pc_fieldsplit_off_diag_use_amat -stokes_ksp_type pipegcr -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type bjacobi -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type bjacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short -stokes_ksp_view

    1912:    test:
    1913:       suffix: 6
    1914:       nsize: 8
    1915:       args: -stokes_ksp_view -stokes_pc_type mg -stokes_pc_mg_levels 2 -stokes_mg_coarse_pc_type telescope -stokes_mg_coarse_pc_telescope_reduction_factor 2 -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_coarse_pc_telescope_subcomm_type contiguous

    1917:    test:
    1918:       suffix: bjacobi
    1919:       nsize: 4
    1920:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason

    1922:    test:
    1923:       suffix: bjacobi_baij
    1924:       nsize: 4
    1925:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
    1926:       output_file: output/ex43_bjacobi.out

    1928:    test:
    1929:       suffix: nested_gmg
    1930:       nsize: 4
    1931:       args: -stokes_pc_fieldsplit_off_diag_use_amat -mx 16 -my 16 -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_fieldsplit_u_pc_type mg -stokes_fieldsplit_u_pc_mg_levels 5 -stokes_fieldsplit_u_pc_mg_galerkin pmat -stokes_fieldsplit_u_ksp_type cg -stokes_fieldsplit_u_ksp_rtol 1.0e-4 -stokes_fieldsplit_u_mg_levels_pc_type jacobi -solcx_eta0 1.0e4 -stokes_fieldsplit_u_ksp_converged_reason -stokes_ksp_converged_reason -stokes_fieldsplit_p_sub_pc_factor_zeropivot 1.e-8

    1933:    test:
    1934:       suffix: fetidp
    1935:       nsize: 8
    1936:       args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_converged_reason -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1

    1938:    test:
    1939:       suffix: fetidp_unsym
    1940:       nsize: 8
    1941:       args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_monitor_true_residual -stokes_ksp_converged_reason -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd

    1943:    test:
    1944:       suffix: bddc_stokes_deluxe
    1945:       nsize: 8
    1946:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -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

    1948:    test:
    1949:       suffix: bddc_stokes_subdomainjump_deluxe
    1950:       nsize: 9
    1951:       args: -c_str 4 -jump_magnitude 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -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

    1953: TEST*/