Actual source code: ex49.c

      1: static char help[] = "   Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
      2:    Material properties E (Youngs modulus) and nu (Poisson ratio) may vary as a function of space. \n\
      3:    The model utilises boundary conditions which produce compression in the x direction. \n\
      4: Options: \n"
      5:                      "\
      6:      -mx : number of elements in x-direction \n\
      7:      -my : number of elements in y-direction \n\
      8:      -c_str : structure of the coefficients to use. \n"
      9:                      "\
     10:           -c_str 0 => isotropic material with constant coefficients. \n\
     11:                          Parameters: \n\
     12:                              -iso_E  : Youngs modulus \n\
     13:                              -iso_nu : Poisson ratio \n\
     14:           -c_str 1 => step function in the material properties in x. \n\
     15:                          Parameters: \n\
     16:                               -step_E0  : Youngs modulus to the left of the step \n\
     17:                               -step_nu0 : Poisson ratio to the left of the step \n\
     18:                               -step_E1  : Youngs modulus to the right of the step \n\
     19:                               -step_n1  : Poisson ratio to the right of the step \n\
     20:                               -step_xc  : x coordinate of the step \n"
     21:                      "\
     22:           -c_str 2 => checkerboard material with alternating properties. \n\
     23:                       Repeats the following pattern throughout the domain. For example with 4 materials specified, we would heve \n\
     24:                       -------------------------\n\
     25:                       |  D  |  A  |  B  |  C  |\n\
     26:                       ------|-----|-----|------\n\
     27:                       |  C  |  D  |  A  |  B  |\n\
     28:                       ------|-----|-----|------\n\
     29:                       |  B  |  C  |  D  |  A  |\n\
     30:                       ------|-----|-----|------\n\
     31:                       |  A  |  B  |  C  |  D  |\n\
     32:                       -------------------------\n\
     33:                       \n\
     34:                          Parameters: \n\
     35:                               -brick_E    : a comma separated list of Young's modulii \n\
     36:                               -brick_nu   : a comma separated list of Poisson ratios  \n\
     37:                               -brick_span : the number of elements in x and y each brick will span \n\
     38:           -c_str 3 => sponge-like material with alternating properties. \n\
     39:                       Repeats the following pattern throughout the domain \n"
     40:                      "\
     41:                       -----------------------------\n\
     42:                       |       [background]        |\n\
     43:                       |          E0,nu0           |\n\
     44:                       |     -----------------     |\n\
     45:                       |     |  [inclusion]  |     |\n\
     46:                       |     |    E1,nu1     |     |\n\
     47:                       |     |               |     |\n\
     48:                       |     | <---- w ----> |     |\n\
     49:                       |     |               |     |\n\
     50:                       |     |               |     |\n\
     51:                       |     -----------------     |\n\
     52:                       |                           |\n\
     53:                       |                           |\n\
     54:                       -----------------------------\n\
     55:                       <--------  t + w + t ------->\n\
     56:                       \n\
     57:                          Parameters: \n\
     58:                               -sponge_E0  : Youngs modulus of the surrounding material \n\
     59:                               -sponge_E1  : Youngs modulus of the inclusion \n\
     60:                               -sponge_nu0 : Poisson ratio of the surrounding material \n\
     61:                               -sponge_nu1 : Poisson ratio of the inclusion \n\
     62:                               -sponge_t   : the number of elements defining the border around each inclusion \n\
     63:                               -sponge_w   : the number of elements in x and y each inclusion will span\n\
     64:      -use_gp_coords : Evaluate the Youngs modulus, Poisson ratio and the body force at the global coordinates of the quadrature points.\n\
     65:      By default, E, nu and the body force are evaluated at the element center and applied as a constant over the entire element.\n\
     66:      -use_nonsymbc : Option to use non-symmetric boundary condition imposition. This choice will use less memory.";

     68: /* Contributed by Dave May */

     70: #include <petscksp.h>
     71: #include <petscdm.h>
     72: #include <petscdmda.h>

     74: static PetscErrorCode DMDABCApplyCompression(DM, Mat, Vec);
     75: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff);

     77: #define NSD          2 /* number of spatial dimensions */
     78: #define NODES_PER_EL 4 /* nodes per element */
     79: #define U_DOFS       2 /* degrees of freedom per displacement node */
     80: #define GAUSS_POINTS 4

     82: /* cell based evaluation */
     83: typedef struct {
     84:   PetscScalar E, nu, fx, fy;
     85: } Coefficients;

     87: /* Gauss point based evaluation 8+4+4+4 = 20 */
     88: typedef struct {
     89:   PetscScalar gp_coords[2 * GAUSS_POINTS];
     90:   PetscScalar E[GAUSS_POINTS];
     91:   PetscScalar nu[GAUSS_POINTS];
     92:   PetscScalar fx[GAUSS_POINTS];
     93:   PetscScalar fy[GAUSS_POINTS];
     94: } GaussPointCoefficients;

     96: typedef struct {
     97:   PetscScalar ux_dof;
     98:   PetscScalar uy_dof;
     99: } ElasticityDOF;

    101: /*

    103:  D = E/((1+nu)(1-2nu)) * [ 1-nu   nu        0     ]
    104:                          [  nu   1-nu       0     ]
    105:                          [  0     0   0.5*(1-2nu) ]

    107:  B = [ d_dx   0   ]
    108:      [  0    d_dy ]
    109:      [ d_dy  d_dx ]

    111:  */

    113: /* FEM routines */
    114: /*
    115:  Element: Local basis function ordering
    116:  1-----2
    117:  |     |
    118:  |     |
    119:  0-----3
    120:  */
    121: static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
    122: {
    123:   PetscScalar xi  = _xi[0];
    124:   PetscScalar eta = _xi[1];

    126:   Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
    127:   Ni[1] = 0.25 * (1.0 - xi) * (1.0 + eta);
    128:   Ni[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
    129:   Ni[3] = 0.25 * (1.0 + xi) * (1.0 - eta);
    130: }

    132: static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
    133: {
    134:   PetscScalar xi  = _xi[0];
    135:   PetscScalar eta = _xi[1];

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

    142:   GNi[1][0] = -0.25 * (1.0 - xi);
    143:   GNi[1][1] = 0.25 * (1.0 - xi);
    144:   GNi[1][2] = 0.25 * (1.0 + xi);
    145:   GNi[1][3] = -0.25 * (1.0 + xi);
    146: }

    148: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
    149: {
    150:   PetscScalar J00, J01, J10, J11, J;
    151:   PetscScalar iJ00, iJ01, iJ10, iJ11;
    152:   PetscInt    i;

    154:   J00 = J01 = J10 = J11 = 0.0;
    155:   for (i = 0; i < NODES_PER_EL; i++) {
    156:     PetscScalar cx = coords[2 * i + 0];
    157:     PetscScalar cy = coords[2 * i + 1];

    159:     J00 = J00 + GNi[0][i] * cx; /* J_xx = dx/dxi */
    160:     J01 = J01 + GNi[0][i] * cy; /* J_xy = dy/dxi */
    161:     J10 = J10 + GNi[1][i] * cx; /* J_yx = dx/deta */
    162:     J11 = J11 + GNi[1][i] * cy; /* J_yy = dy/deta */
    163:   }
    164:   J = (J00 * J11) - (J01 * J10);

    166:   iJ00 = J11 / J;
    167:   iJ01 = -J01 / J;
    168:   iJ10 = -J10 / J;
    169:   iJ11 = J00 / J;

    171:   for (i = 0; i < NODES_PER_EL; i++) {
    172:     GNx[0][i] = GNi[0][i] * iJ00 + GNi[1][i] * iJ01;
    173:     GNx[1][i] = GNi[0][i] * iJ10 + GNi[1][i] * iJ11;
    174:   }

    176:   if (det_J) *det_J = J;
    177: }

    179: static void ConstructGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
    180: {
    181:   *ngp         = 4;
    182:   gp_xi[0][0]  = -0.57735026919;
    183:   gp_xi[0][1]  = -0.57735026919;
    184:   gp_xi[1][0]  = -0.57735026919;
    185:   gp_xi[1][1]  = 0.57735026919;
    186:   gp_xi[2][0]  = 0.57735026919;
    187:   gp_xi[2][1]  = 0.57735026919;
    188:   gp_xi[3][0]  = 0.57735026919;
    189:   gp_xi[3][1]  = -0.57735026919;
    190:   gp_weight[0] = 1.0;
    191:   gp_weight[1] = 1.0;
    192:   gp_weight[2] = 1.0;
    193:   gp_weight[3] = 1.0;
    194: }

    196: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da, PetscInt **_lx, PetscInt **_ly)
    197: {
    198:   PetscMPIInt  rank;
    199:   PetscInt     proc_I, proc_J;
    200:   PetscInt     cpu_x, cpu_y;
    201:   PetscInt     local_mx, local_my;
    202:   Vec          vlx, vly;
    203:   PetscInt    *LX, *LY, i;
    204:   PetscScalar *_a;
    205:   Vec          V_SEQ;
    206:   VecScatter   ctx;

    208:   PetscFunctionBeginUser;
    209:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

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

    213:   proc_J = rank / cpu_x;
    214:   proc_I = rank - cpu_x * proc_J;

    216:   PetscCall(PetscMalloc1(cpu_x, &LX));
    217:   PetscCall(PetscMalloc1(cpu_y, &LY));

    219:   PetscCall(DMDAGetElementsSizes(da, &local_mx, &local_my, NULL));
    220:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vlx));
    221:   PetscCall(VecSetSizes(vlx, PETSC_DECIDE, cpu_x));
    222:   PetscCall(VecSetFromOptions(vlx));

    224:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
    225:   PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
    226:   PetscCall(VecSetFromOptions(vly));

    228:   PetscCall(VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES));
    229:   PetscCall(VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES));
    230:   PetscCall(VecAssemblyBegin(vlx));
    231:   PetscCall(VecAssemblyEnd(vlx));
    232:   PetscCall(VecAssemblyBegin(vly));
    233:   PetscCall(VecAssemblyEnd(vly));

    235:   PetscCall(VecScatterCreateToAll(vlx, &ctx, &V_SEQ));
    236:   PetscCall(VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    237:   PetscCall(VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    238:   PetscCall(VecGetArray(V_SEQ, &_a));
    239:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
    240:   PetscCall(VecRestoreArray(V_SEQ, &_a));
    241:   PetscCall(VecScatterDestroy(&ctx));
    242:   PetscCall(VecDestroy(&V_SEQ));

    244:   PetscCall(VecScatterCreateToAll(vly, &ctx, &V_SEQ));
    245:   PetscCall(VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    246:   PetscCall(VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
    247:   PetscCall(VecGetArray(V_SEQ, &_a));
    248:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
    249:   PetscCall(VecRestoreArray(V_SEQ, &_a));
    250:   PetscCall(VecScatterDestroy(&ctx));
    251:   PetscCall(VecDestroy(&V_SEQ));

    253:   *_lx = LX;
    254:   *_ly = LY;

    256:   PetscCall(VecDestroy(&vlx));
    257:   PetscCall(VecDestroy(&vly));
    258:   PetscFunctionReturn(PETSC_SUCCESS);
    259: }

    261: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
    262: {
    263:   DM           cda;
    264:   Vec          coords;
    265:   DMDACoor2d **_coords;
    266:   PetscInt     si, sj, nx, ny, i, j;
    267:   FILE        *fp;
    268:   char         fname[PETSC_MAX_PATH_LEN];
    269:   PetscMPIInt  rank;

    271:   PetscFunctionBeginUser;
    272:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    273:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
    274:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
    275:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
    276:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank));

    278:   PetscCall(DMGetCoordinateDM(da, &cda));
    279:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    280:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    281:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    282:   for (j = sj; j < sj + ny - 1; j++) {
    283:     for (i = si; i < si + nx - 1; i++) {
    284:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
    285:       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)));
    286:       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)));
    287:       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)));
    288:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
    289:     }
    290:   }
    291:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

    293:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
    294:   PetscFunctionReturn(PETSC_SUCCESS);
    295: }

    297: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
    298: {
    299:   DM           cda;
    300:   Vec          coords, local_fields;
    301:   DMDACoor2d **_coords;
    302:   FILE        *fp;
    303:   char         fname[PETSC_MAX_PATH_LEN];
    304:   const char  *field_name;
    305:   PetscMPIInt  rank;
    306:   PetscInt     si, sj, nx, ny, i, j;
    307:   PetscInt     n_dofs, d;
    308:   PetscScalar *_fields;

    310:   PetscFunctionBeginUser;
    311:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    312:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
    313:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
    314:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

    316:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
    317:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
    318:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
    319:   for (d = 0; d < n_dofs; d++) {
    320:     PetscCall(DMDAGetFieldName(da, d, &field_name));
    321:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
    322:   }
    323:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

    325:   PetscCall(DMGetCoordinateDM(da, &cda));
    326:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    327:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    328:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));

    330:   PetscCall(DMCreateLocalVector(da, &local_fields));
    331:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
    332:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
    333:   PetscCall(VecGetArray(local_fields, &_fields));

    335:   for (j = sj; j < sj + ny; j++) {
    336:     for (i = si; i < si + nx; i++) {
    337:       PetscScalar coord_x, coord_y;
    338:       PetscScalar field_d;

    340:       coord_x = _coords[j][i].x;
    341:       coord_y = _coords[j][i].y;

    343:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
    344:       for (d = 0; d < n_dofs; d++) {
    345:         field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
    346:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d)));
    347:       }
    348:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
    349:     }
    350:   }
    351:   PetscCall(VecRestoreArray(local_fields, &_fields));
    352:   PetscCall(VecDestroy(&local_fields));

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

    356:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
    357:   PetscFunctionReturn(PETSC_SUCCESS);
    358: }

    360: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
    361: {
    362:   DM                       cda;
    363:   Vec                      local_fields;
    364:   FILE                    *fp;
    365:   char                     fname[PETSC_MAX_PATH_LEN];
    366:   const char              *field_name;
    367:   PetscMPIInt              rank;
    368:   PetscInt                 si, sj, nx, ny, i, j, p;
    369:   PetscInt                 n_dofs, d;
    370:   GaussPointCoefficients **_coefficients;

    372:   PetscFunctionBeginUser;
    373:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    374:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
    375:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
    376:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

    378:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
    379:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
    380:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
    381:   for (d = 0; d < n_dofs; d++) {
    382:     PetscCall(DMDAGetFieldName(da, d, &field_name));
    383:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
    384:   }
    385:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

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

    390:   PetscCall(DMCreateLocalVector(da, &local_fields));
    391:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
    392:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
    393:   PetscCall(DMDAVecGetArray(da, local_fields, &_coefficients));

    395:   for (j = sj; j < sj + ny; j++) {
    396:     for (i = si; i < si + nx; i++) {
    397:       PetscScalar coord_x, coord_y;

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

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

    405:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e %1.6e %1.6e\n", (double)PetscRealPart(_coefficients[j][i].E[p]), (double)PetscRealPart(_coefficients[j][i].nu[p]), (double)PetscRealPart(_coefficients[j][i].fx[p]),
    406:                                (double)PetscRealPart(_coefficients[j][i].fy[p])));
    407:       }
    408:     }
    409:   }
    410:   PetscCall(DMDAVecRestoreArray(da, local_fields, &_coefficients));
    411:   PetscCall(VecDestroy(&local_fields));

    413:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
    414:   PetscFunctionReturn(PETSC_SUCCESS);
    415: }

    417: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar E[], PetscScalar nu[])
    418: {
    419:   PetscInt    ngp;
    420:   PetscScalar gp_xi[GAUSS_POINTS][2];
    421:   PetscScalar gp_weight[GAUSS_POINTS];
    422:   PetscInt    p, i, j, k, l;
    423:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    424:   PetscScalar J_p;
    425:   PetscScalar B[3][U_DOFS * NODES_PER_EL];
    426:   PetscScalar prop_E, prop_nu, factor, constit_D[3][3];

    428:   /* define quadrature rule */
    429:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    431:   /* evaluate integral */
    432:   for (p = 0; p < ngp; p++) {
    433:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
    434:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);

    436:     for (i = 0; i < NODES_PER_EL; i++) {
    437:       PetscScalar d_dx_i = GNx_p[0][i];
    438:       PetscScalar d_dy_i = GNx_p[1][i];

    440:       B[0][2 * i]     = d_dx_i;
    441:       B[0][2 * i + 1] = 0.0;
    442:       B[1][2 * i]     = 0.0;
    443:       B[1][2 * i + 1] = d_dy_i;
    444:       B[2][2 * i]     = d_dy_i;
    445:       B[2][2 * i + 1] = d_dx_i;
    446:     }

    448:     /* form D for the quadrature point */
    449:     prop_E          = E[p];
    450:     prop_nu         = nu[p];
    451:     factor          = prop_E / ((1.0 + prop_nu) * (1.0 - 2.0 * prop_nu));
    452:     constit_D[0][0] = 1.0 - prop_nu;
    453:     constit_D[0][1] = prop_nu;
    454:     constit_D[0][2] = 0.0;
    455:     constit_D[1][0] = prop_nu;
    456:     constit_D[1][1] = 1.0 - prop_nu;
    457:     constit_D[1][2] = 0.0;
    458:     constit_D[2][0] = 0.0;
    459:     constit_D[2][1] = 0.0;
    460:     constit_D[2][2] = 0.5 * (1.0 - 2.0 * prop_nu);
    461:     for (i = 0; i < 3; i++) {
    462:       for (j = 0; j < 3; j++) constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
    463:     }

    465:     /* form Bt tildeD B */
    466:     /*
    467:      Ke_ij = Bt_ik . D_kl . B_lj
    468:      = B_ki . D_kl . B_lj
    469:      */
    470:     for (i = 0; i < 8; i++) {
    471:       for (j = 0; j < 8; j++) {
    472:         for (k = 0; k < 3; k++) {
    473:           for (l = 0; l < 3; l++) Ke[8 * i + j] = Ke[8 * i + j] + B[k][i] * constit_D[k][l] * B[l][j];
    474:         }
    475:       }
    476:     }

    478:   } /* end quadrature */
    479: }

    481: static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
    482: {
    483:   PetscInt    ngp;
    484:   PetscScalar gp_xi[GAUSS_POINTS][2];
    485:   PetscScalar gp_weight[GAUSS_POINTS];
    486:   PetscInt    p, i;
    487:   PetscScalar Ni_p[NODES_PER_EL];
    488:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
    489:   PetscScalar J_p, fac;

    491:   /* define quadrature rule */
    492:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

    494:   /* evaluate integral */
    495:   for (p = 0; p < ngp; p++) {
    496:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
    497:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
    498:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
    499:     fac = gp_weight[p] * J_p;

    501:     for (i = 0; i < NODES_PER_EL; i++) {
    502:       Fe[NSD * i] += fac * Ni_p[i] * fx[p];
    503:       Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
    504:     }
    505:   }
    506: }

    508: /*
    509:  i,j are the element indices
    510:  The unknown is a vector quantity.
    511:  The s[].c is used to indicate the degree of freedom.
    512:  */
    513: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[], PetscInt i, PetscInt j)
    514: {
    515:   PetscFunctionBeginUser;
    516:   /* displacement */
    517:   /* node 0 */
    518:   s_u[0].i = i;
    519:   s_u[0].j = j;
    520:   s_u[0].c = 0; /* Ux0 */
    521:   s_u[1].i = i;
    522:   s_u[1].j = j;
    523:   s_u[1].c = 1; /* Uy0 */

    525:   /* node 1 */
    526:   s_u[2].i = i;
    527:   s_u[2].j = j + 1;
    528:   s_u[2].c = 0; /* Ux1 */
    529:   s_u[3].i = i;
    530:   s_u[3].j = j + 1;
    531:   s_u[3].c = 1; /* Uy1 */

    533:   /* node 2 */
    534:   s_u[4].i = i + 1;
    535:   s_u[4].j = j + 1;
    536:   s_u[4].c = 0; /* Ux2 */
    537:   s_u[5].i = i + 1;
    538:   s_u[5].j = j + 1;
    539:   s_u[5].c = 1; /* Uy2 */

    541:   /* node 3 */
    542:   s_u[6].i = i + 1;
    543:   s_u[6].j = j;
    544:   s_u[6].c = 0; /* Ux3 */
    545:   s_u[7].i = i + 1;
    546:   s_u[7].j = j;
    547:   s_u[7].c = 1; /* Uy3 */
    548:   PetscFunctionReturn(PETSC_SUCCESS);
    549: }

    551: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
    552: {
    553:   PetscFunctionBeginUser;
    554:   /* get coords for the element */
    555:   el_coords[NSD * 0 + 0] = _coords[ej][ei].x;
    556:   el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
    557:   el_coords[NSD * 1 + 0] = _coords[ej + 1][ei].x;
    558:   el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
    559:   el_coords[NSD * 2 + 0] = _coords[ej + 1][ei + 1].x;
    560:   el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
    561:   el_coords[NSD * 3 + 0] = _coords[ej][ei + 1].x;
    562:   el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
    563:   PetscFunctionReturn(PETSC_SUCCESS);
    564: }

    566: static PetscErrorCode AssembleA_Elasticity(Mat A, DM elas_da, DM properties_da, Vec properties)
    567: {
    568:   DM                       cda;
    569:   Vec                      coords;
    570:   DMDACoor2d             **_coords;
    571:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    572:   PetscInt                 sex, sey, mx, my;
    573:   PetscInt                 ei, ej;
    574:   PetscScalar              Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
    575:   PetscScalar              el_coords[NODES_PER_EL * NSD];
    576:   Vec                      local_properties;
    577:   GaussPointCoefficients **props;
    578:   PetscScalar             *prop_E, *prop_nu;

    580:   PetscFunctionBeginUser;
    581:   /* setup for coords */
    582:   PetscCall(DMGetCoordinateDM(elas_da, &cda));
    583:   PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
    584:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    586:   /* setup for coefficients */
    587:   PetscCall(DMCreateLocalVector(properties_da, &local_properties));
    588:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
    589:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
    590:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

    592:   PetscCall(DMDAGetElementsCorners(elas_da, &sex, &sey, 0));
    593:   PetscCall(DMDAGetElementsSizes(elas_da, &mx, &my, 0));
    594:   for (ej = sey; ej < sey + my; ej++) {
    595:     for (ei = sex; ei < sex + mx; ei++) {
    596:       /* get coords for the element */
    597:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

    599:       /* get coefficients for the element */
    600:       prop_E  = props[ej][ei].E;
    601:       prop_nu = props[ej][ei].nu;

    603:       /* initialise element stiffness matrix */
    604:       PetscCall(PetscMemzero(Ae, sizeof(Ae)));

    606:       /* form element stiffness matrix */
    607:       FormStressOperatorQ1(Ae, el_coords, prop_E, prop_nu);

    609:       /* insert element matrix into global matrix */
    610:       PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));
    611:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
    612:     }
    613:   }
    614:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    615:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

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

    619:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
    620:   PetscCall(VecDestroy(&local_properties));
    621:   PetscFunctionReturn(PETSC_SUCCESS);
    622: }

    624: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F, MatStencil u_eqn[], PetscScalar Fe_u[])
    625: {
    626:   PetscInt n;

    628:   PetscFunctionBeginUser;
    629:   for (n = 0; n < 4; n++) {
    630:     fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].ux_dof         = fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].ux_dof + Fe_u[2 * n];
    631:     fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].uy_dof = fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].uy_dof + Fe_u[2 * n + 1];
    632:   }
    633:   PetscFunctionReturn(PETSC_SUCCESS);
    634: }

    636: static PetscErrorCode AssembleF_Elasticity(Vec F, DM elas_da, DM properties_da, Vec properties)
    637: {
    638:   DM                       cda;
    639:   Vec                      coords;
    640:   DMDACoor2d             **_coords;
    641:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
    642:   PetscInt                 sex, sey, mx, my;
    643:   PetscInt                 ei, ej;
    644:   PetscScalar              Fe[NODES_PER_EL * U_DOFS];
    645:   PetscScalar              el_coords[NODES_PER_EL * NSD];
    646:   Vec                      local_properties;
    647:   GaussPointCoefficients **props;
    648:   PetscScalar             *prop_fx, *prop_fy;
    649:   Vec                      local_F;
    650:   ElasticityDOF          **ff;

    652:   PetscFunctionBeginUser;
    653:   /* setup for coords */
    654:   PetscCall(DMGetCoordinateDM(elas_da, &cda));
    655:   PetscCall(DMGetCoordinatesLocal(elas_da, &coords));
    656:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

    658:   /* setup for coefficients */
    659:   PetscCall(DMGetLocalVector(properties_da, &local_properties));
    660:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
    661:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
    662:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

    664:   /* get access to the vector */
    665:   PetscCall(DMGetLocalVector(elas_da, &local_F));
    666:   PetscCall(VecZeroEntries(local_F));
    667:   PetscCall(DMDAVecGetArray(elas_da, local_F, &ff));

    669:   PetscCall(DMDAGetElementsCorners(elas_da, &sex, &sey, 0));
    670:   PetscCall(DMDAGetElementsSizes(elas_da, &mx, &my, 0));
    671:   for (ej = sey; ej < sey + my; ej++) {
    672:     for (ei = sex; ei < sex + mx; ei++) {
    673:       /* get coords for the element */
    674:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

    676:       /* get coefficients for the element */
    677:       prop_fx = props[ej][ei].fx;
    678:       prop_fy = props[ej][ei].fy;

    680:       /* initialise element stiffness matrix */
    681:       PetscCall(PetscMemzero(Fe, sizeof(Fe)));

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

    686:       /* insert element matrix into global matrix */
    687:       PetscCall(DMDAGetElementEqnums_u(u_eqn, ei, ej));

    689:       PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, Fe));
    690:     }
    691:   }

    693:   PetscCall(DMDAVecRestoreArray(elas_da, local_F, &ff));
    694:   PetscCall(DMLocalToGlobalBegin(elas_da, local_F, ADD_VALUES, F));
    695:   PetscCall(DMLocalToGlobalEnd(elas_da, local_F, ADD_VALUES, F));
    696:   PetscCall(DMRestoreLocalVector(elas_da, &local_F));

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

    700:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
    701:   PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
    702:   PetscFunctionReturn(PETSC_SUCCESS);
    703: }

    705: static PetscErrorCode solve_elasticity_2d(PetscInt mx, PetscInt my)
    706: {
    707:   DM                       elas_da, da_prop;
    708:   PetscInt                 u_dof, dof, stencil_width;
    709:   Mat                      A;
    710:   PetscInt                 mxl, myl;
    711:   DM                       prop_cda, vel_cda;
    712:   Vec                      prop_coords, vel_coords;
    713:   PetscInt                 si, sj, nx, ny, i, j, p;
    714:   Vec                      f, X;
    715:   PetscInt                 prop_dof, prop_stencil_width;
    716:   Vec                      properties, l_properties;
    717:   MatNullSpace             matnull;
    718:   PetscReal                dx, dy;
    719:   PetscInt                 M, N;
    720:   DMDACoor2d             **_prop_coords, **_vel_coords;
    721:   GaussPointCoefficients **element_props;
    722:   KSP                      ksp_E;
    723:   PetscInt                 coefficient_structure = 0;
    724:   PetscInt                 cpu_x, cpu_y, *lx = NULL, *ly = NULL;
    725:   PetscBool                use_gp_coords = PETSC_FALSE;
    726:   PetscBool                use_nonsymbc  = PETSC_FALSE;
    727:   PetscBool                no_view       = PETSC_FALSE;
    728:   PetscBool                flg;

    730:   PetscFunctionBeginUser;
    731:   /* Generate the da for velocity and pressure */
    732:   /*
    733:    We use Q1 elements for the temperature.
    734:    FEM has a 9-point stencil (BOX) or connectivity pattern
    735:    Num nodes in each direction is mx+1, my+1
    736:    */
    737:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
    738:   dof           = u_dof;
    739:   stencil_width = 1;
    740:   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, &elas_da));

    742:   PetscCall(DMSetMatType(elas_da, MATAIJ));
    743:   PetscCall(DMSetFromOptions(elas_da));
    744:   PetscCall(DMSetUp(elas_da));

    746:   PetscCall(DMDASetFieldName(elas_da, 0, "Ux"));
    747:   PetscCall(DMDASetFieldName(elas_da, 1, "Uy"));

    749:   /* unit box [0,1] x [0,1] */
    750:   PetscCall(DMDASetUniformCoordinates(elas_da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));

    752:   /* Generate element properties, we will assume all material properties are constant over the element */
    753:   /* local number of elements */
    754:   PetscCall(DMDAGetElementsSizes(elas_da, &mxl, &myl, NULL));

    756:   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
    757:   PetscCall(DMDAGetInfo(elas_da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
    758:   PetscCall(DMDAGetElementOwnershipRanges2d(elas_da, &lx, &ly));

    760:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
    761:   prop_stencil_width = 0;
    762:   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));
    763:   PetscCall(DMSetFromOptions(da_prop));
    764:   PetscCall(DMSetUp(da_prop));

    766:   PetscCall(PetscFree(lx));
    767:   PetscCall(PetscFree(ly));

    769:   /* define centroid positions */
    770:   PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    771:   dx = 1.0 / (PetscReal)M;
    772:   dy = 1.0 / (PetscReal)N;

    774:   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, 1.0));

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

    779:   PetscCall(DMCreateGlobalVector(da_prop, &properties));
    780:   PetscCall(DMCreateLocalVector(da_prop, &l_properties));
    781:   PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));

    783:   PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
    784:   PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
    785:   PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));

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

    789:   PetscCall(DMGetCoordinateDM(elas_da, &vel_cda));
    790:   PetscCall(DMGetCoordinatesLocal(elas_da, &vel_coords));
    791:   PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));

    793:   /* interpolate the coordinates */
    794:   for (j = sj; j < sj + ny; j++) {
    795:     for (i = si; i < si + nx; i++) {
    796:       PetscInt    ngp;
    797:       PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
    798:       PetscScalar el_coords[8];

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

    803:       for (p = 0; p < GAUSS_POINTS; p++) {
    804:         PetscScalar gp_x, gp_y;
    805:         PetscInt    n;
    806:         PetscScalar xi_p[2], Ni_p[4];

    808:         xi_p[0] = gp_xi[p][0];
    809:         xi_p[1] = gp_xi[p][1];
    810:         ConstructQ12D_Ni(xi_p, Ni_p);

    812:         gp_x = 0.0;
    813:         gp_y = 0.0;
    814:         for (n = 0; n < NODES_PER_EL; n++) {
    815:           gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
    816:           gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
    817:         }
    818:         element_props[j][i].gp_coords[2 * p]     = gp_x;
    819:         element_props[j][i].gp_coords[2 * p + 1] = gp_y;
    820:       }
    821:     }
    822:   }

    824:   /* define the coefficients */
    825:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, &flg));

    827:   for (j = sj; j < sj + ny; j++) {
    828:     for (i = si; i < si + nx; i++) {
    829:       PetscScalar              centroid_x = _prop_coords[j][i].x; /* centroids of cell */
    830:       PetscScalar              centroid_y = _prop_coords[j][i].y;
    831:       PETSC_UNUSED PetscScalar coord_x, coord_y;

    833:       if (coefficient_structure == 0) { /* isotropic */
    834:         PetscScalar opts_E, opts_nu;

    836:         opts_E  = 1.0;
    837:         opts_nu = 0.33;
    838:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_E", &opts_E, &flg));
    839:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-iso_nu", &opts_nu, &flg));

    841:         for (p = 0; p < GAUSS_POINTS; p++) {
    842:           element_props[j][i].E[p]  = opts_E;
    843:           element_props[j][i].nu[p] = opts_nu;

    845:           element_props[j][i].fx[p] = 0.0;
    846:           element_props[j][i].fy[p] = 0.0;
    847:         }
    848:       } else if (coefficient_structure == 1) { /* step */
    849:         PetscScalar opts_E0, opts_nu0, opts_xc;
    850:         PetscScalar opts_E1, opts_nu1;

    852:         opts_E0 = opts_E1 = 1.0;
    853:         opts_nu0 = opts_nu1 = 0.333;
    854:         opts_xc             = 0.5;
    855:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E0", &opts_E0, &flg));
    856:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu0", &opts_nu0, &flg));
    857:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_E1", &opts_E1, &flg));
    858:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_nu1", &opts_nu1, &flg));
    859:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-step_xc", &opts_xc, &flg));

    861:         for (p = 0; p < GAUSS_POINTS; p++) {
    862:           coord_x = centroid_x;
    863:           coord_y = centroid_y;
    864:           if (use_gp_coords) {
    865:             coord_x = element_props[j][i].gp_coords[2 * p];
    866:             coord_y = element_props[j][i].gp_coords[2 * p + 1];
    867:           }

    869:           element_props[j][i].E[p]  = opts_E0;
    870:           element_props[j][i].nu[p] = opts_nu0;
    871:           if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
    872:             element_props[j][i].E[p]  = opts_E1;
    873:             element_props[j][i].nu[p] = opts_nu1;
    874:           }

    876:           element_props[j][i].fx[p] = 0.0;
    877:           element_props[j][i].fy[p] = 0.0;
    878:         }
    879:       } else if (coefficient_structure == 2) { /* brick */
    880:         PetscReal values_E[10];
    881:         PetscReal values_nu[10];
    882:         PetscInt  nbricks, maxnbricks;
    883:         PetscInt  index, span;
    884:         PetscInt  jj;

    886:         flg        = PETSC_FALSE;
    887:         maxnbricks = 10;
    888:         PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_E", values_E, &maxnbricks, &flg));
    889:         nbricks = maxnbricks;
    890:         PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of E values for each brick");

    892:         flg        = PETSC_FALSE;
    893:         maxnbricks = 10;
    894:         PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-brick_nu", values_nu, &maxnbricks, &flg));
    895:         PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply a list of nu values for each brick");
    896:         PetscCheck(maxnbricks == nbricks, PETSC_COMM_SELF, PETSC_ERR_USER, "User must supply equal numbers of values for E and nu");

    898:         span = 1;
    899:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-brick_span", &span, &flg));

    901:         /* cycle through the indices so that no two material properties are repeated in lines of x or y */
    902:         jj    = (j / span) % nbricks;
    903:         index = (jj + i / span) % nbricks;
    904:         /*printf("j=%d: index = %d \n", j,index); */

    906:         for (p = 0; p < GAUSS_POINTS; p++) {
    907:           element_props[j][i].E[p]  = values_E[index];
    908:           element_props[j][i].nu[p] = values_nu[index];
    909:         }
    910:       } else if (coefficient_structure == 3) { /* sponge */
    911:         PetscScalar opts_E0, opts_nu0;
    912:         PetscScalar opts_E1, opts_nu1;
    913:         PetscInt    opts_t, opts_w;
    914:         PetscInt    ii, jj, ci, cj;

    916:         opts_E0 = opts_E1 = 1.0;
    917:         opts_nu0 = opts_nu1 = 0.333;
    918:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E0", &opts_E0, &flg));
    919:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu0", &opts_nu0, &flg));
    920:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_E1", &opts_E1, &flg));
    921:         PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sponge_nu1", &opts_nu1, &flg));

    923:         opts_t = opts_w = 1;
    924:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_t", &opts_t, &flg));
    925:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-sponge_w", &opts_w, &flg));

    927:         ii = (i) / (opts_t + opts_w + opts_t);
    928:         jj = (j) / (opts_t + opts_w + opts_t);

    930:         ci = i - ii * (opts_t + opts_w + opts_t);
    931:         cj = j - jj * (opts_t + opts_w + opts_t);

    933:         for (p = 0; p < GAUSS_POINTS; p++) {
    934:           element_props[j][i].E[p]  = opts_E0;
    935:           element_props[j][i].nu[p] = opts_nu0;
    936:         }
    937:         if ((ci >= opts_t) && (ci < opts_t + opts_w)) {
    938:           if ((cj >= opts_t) && (cj < opts_t + opts_w)) {
    939:             for (p = 0; p < GAUSS_POINTS; p++) {
    940:               element_props[j][i].E[p]  = opts_E1;
    941:               element_props[j][i].nu[p] = opts_nu1;
    942:             }
    943:           }
    944:         }
    945:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
    946:     }
    947:   }
    948:   PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));

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

    952:   PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
    953:   PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
    954:   PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));

    956:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
    957:   if (!no_view) {
    958:     PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for elasticity eqn.", "properties"));
    959:     PetscCall(DMDACoordViewGnuplot2d(elas_da, "mesh"));
    960:   }

    962:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
    963:   PetscCall(DMCreateMatrix(elas_da, &A));
    964:   PetscCall(DMGetCoordinates(elas_da, &vel_coords));
    965:   PetscCall(MatNullSpaceCreateRigidBody(vel_coords, &matnull));
    966:   PetscCall(MatSetNearNullSpace(A, matnull));
    967:   PetscCall(MatNullSpaceDestroy(&matnull));
    968:   PetscCall(MatCreateVecs(A, &f, &X));

    970:   /* assemble A11 */
    971:   PetscCall(MatZeroEntries(A));
    972:   PetscCall(VecZeroEntries(f));

    974:   PetscCall(AssembleA_Elasticity(A, elas_da, da_prop, properties));
    975:   /* build force vector */
    976:   PetscCall(AssembleF_Elasticity(f, elas_da, da_prop, properties));

    978:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_E));
    979:   PetscCall(KSPSetOptionsPrefix(ksp_E, "elas_")); /* elasticity */

    981:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_nonsymbc", &use_nonsymbc, &flg));
    982:   /* solve */
    983:   if (!use_nonsymbc) {
    984:     Mat        AA;
    985:     Vec        ff, XX;
    986:     IS         is;
    987:     VecScatter scat;

    989:     PetscCall(DMDABCApplySymmetricCompression(elas_da, A, f, &is, &AA, &ff));
    990:     PetscCall(VecDuplicate(ff, &XX));

    992:     PetscCall(KSPSetOperators(ksp_E, AA, AA));
    993:     PetscCall(KSPSetFromOptions(ksp_E));

    995:     PetscCall(KSPSolve(ksp_E, ff, XX));

    997:     /* push XX back into X */
    998:     PetscCall(DMDABCApplyCompression(elas_da, NULL, X));

    1000:     PetscCall(VecScatterCreate(XX, NULL, X, is, &scat));
    1001:     PetscCall(VecScatterBegin(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
    1002:     PetscCall(VecScatterEnd(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD));
    1003:     PetscCall(VecScatterDestroy(&scat));

    1005:     PetscCall(MatDestroy(&AA));
    1006:     PetscCall(VecDestroy(&ff));
    1007:     PetscCall(VecDestroy(&XX));
    1008:     PetscCall(ISDestroy(&is));
    1009:   } else {
    1010:     PetscCall(DMDABCApplyCompression(elas_da, A, f));

    1012:     PetscCall(KSPSetOperators(ksp_E, A, A));
    1013:     PetscCall(KSPSetFromOptions(ksp_E));

    1015:     PetscCall(KSPSolve(ksp_E, f, X));
    1016:   }

    1018:   if (!no_view) PetscCall(DMDAViewGnuplot2d(elas_da, X, "Displacement solution for elasticity eqn.", "X"));
    1019:   PetscCall(KSPDestroy(&ksp_E));

    1021:   PetscCall(VecDestroy(&X));
    1022:   PetscCall(VecDestroy(&f));
    1023:   PetscCall(MatDestroy(&A));

    1025:   PetscCall(DMDestroy(&elas_da));
    1026:   PetscCall(DMDestroy(&da_prop));

    1028:   PetscCall(VecDestroy(&properties));
    1029:   PetscCall(VecDestroy(&l_properties));
    1030:   PetscFunctionReturn(PETSC_SUCCESS);
    1031: }

    1033: int main(int argc, char **args)
    1034: {
    1035:   PetscInt mx, my;

    1037:   PetscFunctionBeginUser;
    1038:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
    1039:   mx = my = 10;
    1040:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
    1041:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
    1042:   PetscCall(solve_elasticity_2d(mx, my));
    1043:   PetscCall(PetscFinalize());
    1044:   return 0;
    1045: }

    1047: /* -------------------------- helpers for boundary conditions -------------------------------- */

    1049: static PetscErrorCode BCApply_EAST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
    1050: {
    1051:   DM                     cda;
    1052:   Vec                    coords;
    1053:   PetscInt               si, sj, nx, ny, i, j;
    1054:   PetscInt               M, N;
    1055:   DMDACoor2d           **_coords;
    1056:   const PetscInt        *g_idx;
    1057:   PetscInt              *bc_global_ids;
    1058:   PetscScalar           *bc_vals;
    1059:   PetscInt               nbcs;
    1060:   PetscInt               n_dofs;
    1061:   ISLocalToGlobalMapping ltogm;

    1063:   PetscFunctionBeginUser;
    1064:   /* enforce bc's */
    1065:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1066:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1068:   PetscCall(DMGetCoordinateDM(da, &cda));
    1069:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1070:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1071:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1072:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1074:   /* --- */

    1076:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
    1077:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

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

    1082:   i = nx - 1;
    1083:   for (j = 0; j < ny; j++) {
    1084:     PetscInt                 local_id;
    1085:     PETSC_UNUSED PetscScalar coordx, coordy;

    1087:     local_id = i + j * nx;

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

    1091:     coordx = _coords[j + sj][i + si].x;
    1092:     coordy = _coords[j + sj][i + si].y;

    1094:     bc_vals[j] = bc_val;
    1095:   }
    1096:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1097:   nbcs = 0;
    1098:   if ((si + nx) == (M)) nbcs = ny;

    1100:   if (b) {
    1101:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1102:     PetscCall(VecAssemblyBegin(b));
    1103:     PetscCall(VecAssemblyEnd(b));
    1104:   }
    1105:   if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1107:   PetscCall(PetscFree(bc_vals));
    1108:   PetscCall(PetscFree(bc_global_ids));

    1110:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1111:   PetscFunctionReturn(PETSC_SUCCESS);
    1112: }

    1114: static PetscErrorCode BCApply_WEST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
    1115: {
    1116:   DM                     cda;
    1117:   Vec                    coords;
    1118:   PetscInt               si, sj, nx, ny, i, j;
    1119:   PetscInt               M, N;
    1120:   DMDACoor2d           **_coords;
    1121:   const PetscInt        *g_idx;
    1122:   PetscInt              *bc_global_ids;
    1123:   PetscScalar           *bc_vals;
    1124:   PetscInt               nbcs;
    1125:   PetscInt               n_dofs;
    1126:   ISLocalToGlobalMapping ltogm;

    1128:   PetscFunctionBeginUser;
    1129:   /* enforce bc's */
    1130:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
    1131:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

    1133:   PetscCall(DMGetCoordinateDM(da, &cda));
    1134:   PetscCall(DMGetCoordinatesLocal(da, &coords));
    1135:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
    1136:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
    1137:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

    1139:   /* --- */

    1141:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
    1142:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

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

    1147:   i = 0;
    1148:   for (j = 0; j < ny; j++) {
    1149:     PetscInt                 local_id;
    1150:     PETSC_UNUSED PetscScalar coordx, coordy;

    1152:     local_id = i + j * nx;

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

    1156:     coordx = _coords[j + sj][i + si].x;
    1157:     coordy = _coords[j + sj][i + si].y;

    1159:     bc_vals[j] = bc_val;
    1160:   }
    1161:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
    1162:   nbcs = 0;
    1163:   if (si == 0) nbcs = ny;

    1165:   if (b) {
    1166:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
    1167:     PetscCall(VecAssemblyBegin(b));
    1168:     PetscCall(VecAssemblyEnd(b));
    1169:   }
    1170:   if (A) PetscCall(MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0));

    1172:   PetscCall(PetscFree(bc_vals));
    1173:   PetscCall(PetscFree(bc_global_ids));

    1175:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
    1176:   PetscFunctionReturn(PETSC_SUCCESS);
    1177: }

    1179: static PetscErrorCode DMDABCApplyCompression(DM elas_da, Mat A, Vec f)
    1180: {
    1181:   PetscFunctionBeginUser;
    1182:   PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, f));
    1183:   PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, f));
    1184:   PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, f));
    1185:   PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, f));
    1186:   PetscFunctionReturn(PETSC_SUCCESS);
    1187: }

    1189: static PetscErrorCode Orthogonalize(PetscInt n, Vec *vecs)
    1190: {
    1191:   PetscInt    i, j;
    1192:   PetscScalar dot;

    1194:   PetscFunctionBegin;
    1195:   for (i = 0; i < n; i++) {
    1196:     PetscCall(VecNormalize(vecs[i], NULL));
    1197:     for (j = i + 1; j < n; j++) {
    1198:       PetscCall(VecDot(vecs[i], vecs[j], &dot));
    1199:       PetscCall(VecAXPY(vecs[j], -dot, vecs[i]));
    1200:     }
    1201:   }
    1202:   PetscFunctionReturn(PETSC_SUCCESS);
    1203: }

    1205: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da, Mat A, Vec f, IS *dofs, Mat *AA, Vec *ff)
    1206: {
    1207:   PetscInt     start, end, m;
    1208:   PetscInt    *unconstrained;
    1209:   PetscInt     cnt, i;
    1210:   Vec          x;
    1211:   PetscScalar *_x;
    1212:   IS           is;
    1213:   VecScatter   scat;

    1215:   PetscFunctionBeginUser;
    1216:   /* push bc's into f and A */
    1217:   PetscCall(VecDuplicate(f, &x));
    1218:   PetscCall(BCApply_EAST(elas_da, 0, -1.0, A, x));
    1219:   PetscCall(BCApply_EAST(elas_da, 1, 0.0, A, x));
    1220:   PetscCall(BCApply_WEST(elas_da, 0, 1.0, A, x));
    1221:   PetscCall(BCApply_WEST(elas_da, 1, 0.0, A, x));

    1223:   /* define which dofs are not constrained */
    1224:   PetscCall(VecGetLocalSize(x, &m));
    1225:   PetscCall(PetscMalloc1(m, &unconstrained));
    1226:   PetscCall(VecGetOwnershipRange(x, &start, &end));
    1227:   PetscCall(VecGetArray(x, &_x));
    1228:   cnt = 0;
    1229:   for (i = 0; i < m; i += 2) {
    1230:     PetscReal val1, val2;

    1232:     val1 = PetscRealPart(_x[i]);
    1233:     val2 = PetscRealPart(_x[i + 1]);
    1234:     if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
    1235:       unconstrained[cnt] = start + i;
    1236:       cnt++;
    1237:       unconstrained[cnt] = start + i + 1;
    1238:       cnt++;
    1239:     }
    1240:   }
    1241:   PetscCall(VecRestoreArray(x, &_x));

    1243:   PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, cnt, unconstrained, PETSC_COPY_VALUES, &is));
    1244:   PetscCall(PetscFree(unconstrained));
    1245:   PetscCall(ISSetBlockSize(is, 2));

    1247:   /* define correction for dirichlet in the rhs */
    1248:   PetscCall(MatMult(A, x, f));
    1249:   PetscCall(VecScale(f, -1.0));

    1251:   /* get new matrix */
    1252:   PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, AA));
    1253:   /* get new vector */
    1254:   PetscCall(MatCreateVecs(*AA, NULL, ff));

    1256:   PetscCall(VecScatterCreate(f, is, *ff, NULL, &scat));
    1257:   PetscCall(VecScatterBegin(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));
    1258:   PetscCall(VecScatterEnd(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD));

    1260:   { /* Constrain near-null space */
    1261:     PetscInt     nvecs;
    1262:     const Vec   *vecs;
    1263:     Vec         *uvecs;
    1264:     PetscBool    has_const;
    1265:     MatNullSpace mnull, unull;

    1267:     PetscCall(MatGetNearNullSpace(A, &mnull));
    1268:     PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvecs, &vecs));
    1269:     PetscCall(VecDuplicateVecs(*ff, nvecs, &uvecs));
    1270:     for (i = 0; i < nvecs; i++) {
    1271:       PetscCall(VecScatterBegin(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
    1272:       PetscCall(VecScatterEnd(scat, vecs[i], uvecs[i], INSERT_VALUES, SCATTER_FORWARD));
    1273:     }
    1274:     PetscCall(Orthogonalize(nvecs, uvecs));
    1275:     PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)A), PETSC_FALSE, nvecs, uvecs, &unull));
    1276:     PetscCall(MatSetNearNullSpace(*AA, unull));
    1277:     PetscCall(MatNullSpaceDestroy(&unull));
    1278:     PetscCall(VecDestroyVecs(nvecs, &uvecs));
    1279:   }

    1281:   PetscCall(VecScatterDestroy(&scat));

    1283:   *dofs = is;
    1284:   PetscCall(VecDestroy(&x));
    1285:   PetscFunctionReturn(PETSC_SUCCESS);
    1286: }

    1288: /*TEST

    1290:    build:
    1291:       requires: !complex !single

    1293:    test:
    1294:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_rtol 5e-3 -elas_ksp_view
    1295:       output_file: output/ex49_1.out

    1297:    test:
    1298:       suffix: 2
    1299:       nsize: 4
    1300:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type gcr -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3

    1302:    test:
    1303:       suffix: 3
    1304:       nsize: 4
    1305:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,10,1000,100 -brick_nu 0.4,0.2,0.3,0.1 -brick_span 3 -elas_pc_type asm -elas_sub_pc_type lu -elas_ksp_rtol 5e-3

    1307:    test:
    1308:       suffix: 4
    1309:       nsize: 4
    1310:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type unpreconditioned -mx 40 -my 40 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_mg_levels_ksp_type chebyshev -elas_pc_type ml -elas_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -elas_mg_levels_pc_type pbjacobi -elas_mg_levels_ksp_max_it 3 -use_nonsymbc -elas_pc_ml_nullspace user
    1311:       requires: ml

    1313:    test:
    1314:       suffix: 5
    1315:       nsize: 3
    1316:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type gamg -elas_mg_fine_ksp_type richardson -elas_mg_fine_pc_type jacobi -elas_mg_fine_pc_jacobi_type rowl1 -elas_mg_fine_pc_jacobi_rowl1_scale .25 -elas_mg_levels_ksp_type chebyshev -elas_mg_levels_ksp_max_it 1 -elas_mg_levels_ksp_chebyshev_esteig 0.2,1.1 -elas_mg_levels_pc_type jacobi -elas_pc_gamg_esteig_ksp_type cg

    1318:    test:
    1319:       suffix: 6
    1320:       nsize: 4
    1321:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type lu

    1323:    test:
    1324:       suffix: 7
    1325:       nsize: 4
    1326:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipegcr -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15

    1328:    test:
    1329:       suffix: 8
    1330:       nsize: 4
    1331:       args: -mx 20 -my 30 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type pipefgmres -elas_pc_type asm -elas_sub_pc_type ksp -elas_sub_ksp_ksp_type cg -elas_sub_ksp_ksp_max_it 15

    1333:    test:
    1334:       suffix: hypre_nullspace
    1335:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
    1336:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type hypre -elas_pc_hypre_boomeramg_nodal_coarsen 6 -elas_pc_hypre_boomeramg_vec_interp_variant 3 -elas_pc_hypre_boomeramg_interp_type ext+i -elas_ksp_view

    1338:    test:
    1339:       nsize: 4
    1340:       suffix: bddc
    1341:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic

    1343:    test:
    1344:       nsize: 4
    1345:       suffix: bddc_unsym
    1346:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0

    1348:    test:
    1349:       nsize: 4
    1350:       suffix: bddc_unsym_deluxe
    1351:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic -use_nonsymbc -elas_pc_bddc_symmetric 0 -elas_pc_bddc_use_deluxe_scaling -elas_sub_schurs_symmetric 0

    1353:    test:
    1354:       nsize: 4
    1355:       suffix: fetidp_unsym_deluxe
    1356:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type fetidp -elas_fetidp_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_fetidp_bddc_pc_bddc_monolithic -use_nonsymbc -elas_fetidp_bddc_pc_bddc_use_deluxe_scaling -elas_fetidp_bddc_sub_schurs_symmetric 0 -elas_fetidp_bddc_pc_bddc_deluxe_singlemat

    1358:    test:
    1359:       nsize: 4
    1360:       suffix: bddc_layerjump
    1361:       args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 3 -sponge_E0 1 -sponge_E1 1000 -sponge_nu0 0.4 -sponge_nu1 0.2 -sponge_t 1 -sponge_w 8 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_ksp_norm_type natural

    1363:    test:
    1364:       nsize: 4
    1365:       suffix: bddc_subdomainjump
    1366:       args: -mx 40 -my 40 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 20 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_is_use_stiffness_scaling -elas_ksp_norm_type natural

    1368:    test:
    1369:       nsize: 9
    1370:       suffix: bddc_subdomainjump_deluxe
    1371:       args: -mx 30 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 10 -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_bddc_use_deluxe_scaling -elas_ksp_norm_type natural -elas_pc_bddc_schur_layers 1
    1372: TEST*/