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 evaulated 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);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    178: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
    179: {
    180:   *ngp         = 4;
    181:   gp_xi[0][0]  = -0.57735026919; gp_xi[0][1] = -0.57735026919;
    182:   gp_xi[1][0]  = -0.57735026919; gp_xi[1][1] =  0.57735026919;
    183:   gp_xi[2][0]  =  0.57735026919; gp_xi[2][1] =  0.57735026919;
    184:   gp_xi[3][0]  =  0.57735026919; 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: {
    199:   /* velocity */
    200:   /* node 0 */
    201:   s_u[0].i = i; s_u[0].j = j; s_u[0].c = 0;       /* Vx0 */
    202:   s_u[1].i = i; s_u[1].j = j; s_u[1].c = 1;       /* Vy0 */

    204:   /* node 1 */
    205:   s_u[2].i = i; s_u[2].j = j+1; s_u[2].c = 0;     /* Vx1 */
    206:   s_u[3].i = i; s_u[3].j = j+1; s_u[3].c = 1;     /* Vy1 */

    208:   /* node 2 */
    209:   s_u[4].i = i+1; s_u[4].j = j+1; s_u[4].c = 0;   /* Vx2 */
    210:   s_u[5].i = i+1; s_u[5].j = j+1; s_u[5].c = 1;   /* Vy2 */

    212:   /* node 3 */
    213:   s_u[6].i = i+1; s_u[6].j = j; s_u[6].c = 0;     /* Vx3 */
    214:   s_u[7].i = i+1; s_u[7].j = j; s_u[7].c = 1;     /* Vy3 */

    216:   /* pressure */
    217:   s_p[0].i = i;   s_p[0].j = j;   s_p[0].c = 2; /* P0 */
    218:   s_p[1].i = i;   s_p[1].j = j+1; s_p[1].c = 2; /* P1 */
    219:   s_p[2].i = i+1; s_p[2].j = j+1; s_p[2].c = 2; /* P2 */
    220:   s_p[3].i = i+1; s_p[3].j = j;   s_p[3].c = 2; /* P3 */
    221:   return(0);
    222: }

    224: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
    225: {
    227:   PetscMPIInt    rank;
    228:   PetscInt       proc_I,proc_J;
    229:   PetscInt       cpu_x,cpu_y;
    230:   PetscInt       local_mx,local_my;
    231:   Vec            vlx,vly;
    232:   PetscInt       *LX,*LY,i;
    233:   PetscScalar    *_a;
    234:   Vec            V_SEQ;
    235:   VecScatter     ctx;

    238:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    240:   DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);

    242:   proc_J = rank/cpu_x;
    243:   proc_I = rank-cpu_x*proc_J;

    245:   PetscMalloc1(cpu_x,&LX);
    246:   PetscMalloc1(cpu_y,&LY);

    248:   DMDAGetElementsSizes(da,&local_mx,&local_my,NULL);
    249:   VecCreate(PETSC_COMM_WORLD,&vlx);
    250:   VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
    251:   VecSetFromOptions(vlx);

    253:   VecCreate(PETSC_COMM_WORLD,&vly);
    254:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
    255:   VecSetFromOptions(vly);

    257:   VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
    258:   VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
    259:   VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
    260:   VecAssemblyBegin(vly);VecAssemblyEnd(vly);

    262:   VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
    263:   VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
    264:   VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
    265:   VecGetArray(V_SEQ,&_a);
    266:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
    267:   VecRestoreArray(V_SEQ,&_a);
    268:   VecScatterDestroy(&ctx);
    269:   VecDestroy(&V_SEQ);

    271:   VecScatterCreateToAll(vly,&ctx,&V_SEQ);
    272:   VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
    273:   VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
    274:   VecGetArray(V_SEQ,&_a);
    275:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
    276:   VecRestoreArray(V_SEQ,&_a);
    277:   VecScatterDestroy(&ctx);
    278:   VecDestroy(&V_SEQ);

    280:   *_lx = LX;
    281:   *_ly = LY;

    283:   VecDestroy(&vlx);
    284:   VecDestroy(&vly);
    285:   return(0);
    286: }

    288: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
    289: {
    290:   DM             cda;
    291:   Vec            coords;
    292:   DMDACoor2d     **_coords;
    293:   PetscInt       si,sj,nx,ny,i,j;
    294:   FILE           *fp;
    295:   char           fname[PETSC_MAX_PATH_LEN];
    296:   PetscMPIInt    rank;

    300:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    301:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
    302:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
    303:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

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

    307:   DMGetCoordinateDM(da,&cda);
    308:   DMGetCoordinatesLocal(da,&coords);
    309:   DMDAVecGetArrayRead(cda,coords,&_coords);
    310:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    311:   for (j = sj; j < sj+ny-1; j++) {
    312:     for (i = si; i < si+nx-1; i++) {
    313:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
    314:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
    315:       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));
    316:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
    317:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
    318:     }
    319:   }
    320:   DMDAVecRestoreArrayRead(cda,coords,&_coords);

    322:   PetscFClose(PETSC_COMM_SELF,fp);
    323:   return(0);
    324: }

    326: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
    327: {
    328:   DM             cda;
    329:   Vec            coords,local_fields;
    330:   DMDACoor2d     **_coords;
    331:   FILE           *fp;
    332:   char           fname[PETSC_MAX_PATH_LEN];
    333:   PetscMPIInt    rank;
    334:   PetscInt       si,sj,nx,ny,i,j;
    335:   PetscInt       n_dofs,d;
    336:   PetscScalar    *_fields;

    340:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    341:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
    342:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
    343:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

    345:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
    346:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
    347:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
    348:   for (d = 0; d < n_dofs; d++) {
    349:     const char *field_name;
    350:     DMDAGetFieldName(da,d,&field_name);
    351:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
    352:   }
    353:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");

    355:   DMGetCoordinateDM(da,&cda);
    356:   DMGetCoordinatesLocal(da,&coords);
    357:   DMDAVecGetArray(cda,coords,&_coords);
    358:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

    360:   DMCreateLocalVector(da,&local_fields);
    361:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
    362:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
    363:   VecGetArray(local_fields,&_fields);

    365:   for (j = sj; j < sj+ny; j++) {
    366:     for (i = si; i < si+nx; i++) {
    367:       PetscScalar coord_x,coord_y;
    368:       PetscScalar field_d;

    370:       coord_x = _coords[j][i].x;
    371:       coord_y = _coords[j][i].y;

    373:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
    374:       for (d = 0; d < n_dofs; d++) {
    375:         field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
    376:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
    377:       }
    378:       PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
    379:     }
    380:   }
    381:   VecRestoreArray(local_fields,&_fields);
    382:   VecDestroy(&local_fields);

    384:   DMDAVecRestoreArray(cda,coords,&_coords);

    386:   PetscFClose(PETSC_COMM_SELF,fp);
    387:   return(0);
    388: }

    390: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
    391: {
    392:   DM                     cda;
    393:   Vec                    local_fields;
    394:   FILE                   *fp;
    395:   char                   fname[PETSC_MAX_PATH_LEN];
    396:   PetscMPIInt            rank;
    397:   PetscInt               si,sj,nx,ny,i,j,p;
    398:   PetscInt               n_dofs,d;
    399:   GaussPointCoefficients **_coefficients;
    400:   PetscErrorCode         ierr;

    403:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    404:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
    405:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
    406:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

    408:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
    409:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
    410:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
    411:   for (d = 0; d < n_dofs; d++) {
    412:     const char *field_name;
    413:     DMDAGetFieldName(da,d,&field_name);
    414:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
    415:   }
    416:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");

    418:   DMGetCoordinateDM(da,&cda);
    419:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

    421:   DMCreateLocalVector(da,&local_fields);
    422:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
    423:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
    424:   DMDAVecGetArray(da,local_fields,&_coefficients);

    426:   for (j = sj; j < sj+ny; j++) {
    427:     for (i = si; i < si+nx; i++) {
    428:       PetscScalar coord_x,coord_y;

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

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

    436:         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]));
    437:         PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
    438:       }
    439:     }
    440:   }
    441:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
    442:   VecDestroy(&local_fields);

    444:   PetscFClose(PETSC_COMM_SELF,fp);
    445:   return(0);
    446: }

    448: 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)
    449: {
    450:   PetscInt ij;
    451:   PetscInt r,c,nc;

    453:   nc = u_NPE*u_dof;
    454:   r = w_dof*wi+wd;
    455:   c = u_dof*ui+ud;
    456:   ij = r*nc+c;
    457:   return ij;
    458: }

    460: /*
    461:  D = [ 2.eta   0   0   ]
    462:      [   0   2.eta 0   ]
    463:      [   0     0   eta ]

    465:  B = [ d_dx   0   ]
    466:      [  0    d_dy ]
    467:      [ d_dy  d_dx ]
    468: */
    469: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
    470: {
    471:   PetscInt       ngp;
    472:   PetscScalar    gp_xi[GAUSS_POINTS][2];
    473:   PetscScalar    gp_weight[GAUSS_POINTS];
    474:   PetscInt       p,i,j,k;
    475:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    476:   PetscScalar    J_p,tildeD[3];
    477:   PetscScalar    B[3][U_DOFS*NODES_PER_EL];

    479:   /* define quadrature rule */
    480:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    482:   /* evaluate integral */
    483:   for (p = 0; p < ngp; p++) {
    484:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
    485:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

    487:     for (i = 0; i < NODES_PER_EL; i++) {
    488:       PetscScalar d_dx_i = GNx_p[0][i];
    489:       PetscScalar d_dy_i = GNx_p[1][i];

    491:       B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
    492:       B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
    493:       B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
    494:     }

    496:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
    497:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
    498:     tildeD[2] =       gp_weight[p]*J_p*eta[p];

    500:     /* form Bt tildeD B */
    501:     /*
    502:     Ke_ij = Bt_ik . D_kl . B_lj
    503:           = B_ki . D_kl . B_lj
    504:           = B_ki . D_kk . B_kj
    505:     */
    506:     for (i = 0; i < 8; i++) {
    507:       for (j = 0; j < 8; j++) {
    508:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
    509:           Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
    510:         }
    511:       }
    512:     }
    513:   }
    514: }

    516: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
    517: {
    518:   PetscInt    ngp;
    519:   PetscScalar gp_xi[GAUSS_POINTS][2];
    520:   PetscScalar gp_weight[GAUSS_POINTS];
    521:   PetscInt    p,i,j,di;
    522:   PetscScalar Ni_p[NODES_PER_EL];
    523:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    524:   PetscScalar J_p,fac;

    526:   /* define quadrature rule */
    527:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    529:   /* evaluate integral */
    530:   for (p = 0; p < ngp; p++) {
    531:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
    532:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
    533:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
    534:     fac = gp_weight[p]*J_p;

    536:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
    537:       for (di = 0; di < NSD; di++) { /* u dofs */
    538:         for (j = 0; j < 4; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
    539:           PetscInt IJ;
    540:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);

    542:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
    543:         }
    544:       }
    545:     }
    546:   }
    547: }

    549: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
    550: {
    551:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
    552:   PetscInt    i,j;
    553:   PetscInt    nr_g,nc_g;

    555:   PetscMemzero(Ge,sizeof(Ge));
    556:   FormGradientOperatorQ1(Ge,coords);

    558:   nr_g = U_DOFS*NODES_PER_EL;
    559:   nc_g = P_DOFS*NODES_PER_EL;

    561:   for (i = 0; i < nr_g; i++) {
    562:     for (j = 0; j < nc_g; j++) {
    563:       De[nr_g*j+i] = Ge[nc_g*i+j];
    564:     }
    565:   }
    566: }

    568: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
    569: {
    570:   PetscInt    ngp;
    571:   PetscScalar gp_xi[GAUSS_POINTS][2];
    572:   PetscScalar gp_weight[GAUSS_POINTS];
    573:   PetscInt    p,i,j;
    574:   PetscScalar Ni_p[NODES_PER_EL];
    575:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    576:   PetscScalar J_p,fac,eta_avg;

    578:   /* define quadrature rule */
    579:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    581:   /* evaluate integral */
    582:   for (p = 0; p < ngp; p++) {
    583:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
    584:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
    585:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
    586:     fac = gp_weight[p]*J_p;

    588:     for (i = 0; i < NODES_PER_EL; i++) {
    589:       for (j = 0; j < NODES_PER_EL; j++) {
    590:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
    591:       }
    592:     }
    593:   }

    595:   /* scale */
    596:   eta_avg = 0.0;
    597:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    598:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
    599:   fac     = 1.0/eta_avg;
    600:   for (i = 0; i < NODES_PER_EL; i++) {
    601:     for (j = 0; j < NODES_PER_EL; j++) {
    602:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
    603:     }
    604:   }
    605: }

    607: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
    608: {
    609:   PetscInt    ngp;
    610:   PetscScalar gp_xi[GAUSS_POINTS][2];
    611:   PetscScalar gp_weight[GAUSS_POINTS];
    612:   PetscInt    p,i,j;
    613:   PetscScalar Ni_p[NODES_PER_EL];
    614:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    615:   PetscScalar J_p,fac,eta_avg;

    617:   /* define quadrature rule */
    618:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    620:   /* evaluate integral */
    621:   for (p = 0; p < ngp; p++) {
    622:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
    623:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
    624:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
    625:     fac = gp_weight[p]*J_p;

    627:     for (i = 0; i < NODES_PER_EL; i++) {
    628:       for (j = 0; j < NODES_PER_EL; j++) {
    629:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
    630:       }
    631:     }
    632:   }

    634:   /* scale */
    635:   eta_avg = 0.0;
    636:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    637:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
    638:   fac     = 1.0/eta_avg;
    639:   for (i = 0; i < NODES_PER_EL; i++) {
    640:     for (j = 0; j < NODES_PER_EL; j++) {
    641:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
    642:     }
    643:   }
    644: }

    646: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
    647: {
    648:   PetscInt    ngp;
    649:   PetscScalar gp_xi[GAUSS_POINTS][2];
    650:   PetscScalar gp_weight[GAUSS_POINTS];
    651:   PetscInt    p,i;
    652:   PetscScalar Ni_p[NODES_PER_EL];
    653:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    654:   PetscScalar J_p,fac;

    656:   /* define quadrature rule */
    657:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    659:   /* evaluate integral */
    660:   for (p = 0; p < ngp; p++) {
    661:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
    662:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
    663:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
    664:     fac = gp_weight[p]*J_p;

    666:     for (i = 0; i < NODES_PER_EL; i++) {
    667:       Fe[NSD*i]   += fac*Ni_p[i]*fx[p];
    668:       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
    669:     }
    670:   }
    671: }

    673: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
    674: {
    676:   /* get coords for the element */
    677:   el_coords[NSD*0] = _coords[ej][ei].x;     el_coords[NSD*0+1] = _coords[ej][ei].y;
    678:   el_coords[NSD*1] = _coords[ej+1][ei].x;   el_coords[NSD*1+1] = _coords[ej+1][ei].y;
    679:   el_coords[NSD*2] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
    680:   el_coords[NSD*3] = _coords[ej][ei+1].x;   el_coords[NSD*3+1] = _coords[ej][ei+1].y;
    681:   return(0);
    682: }

    684: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
    685: {
    686:   DM                     cda;
    687:   Vec                    coords;
    688:   DMDACoor2d             **_coords;
    689:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    690:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
    691:   PetscInt               sex,sey,mx,my;
    692:   PetscInt               ei,ej;
    693:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
    694:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
    695:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
    696:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
    697:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    698:   Vec                    local_properties;
    699:   GaussPointCoefficients **props;
    700:   PetscScalar            *prop_eta;
    701:   PetscErrorCode         ierr;

    704:   /* setup for coords */
    705:   DMGetCoordinateDM(stokes_da,&cda);
    706:   DMGetCoordinatesLocal(stokes_da,&coords);
    707:   DMDAVecGetArray(cda,coords,&_coords);

    709:   /* setup for coefficients */
    710:   DMCreateLocalVector(properties_da,&local_properties);
    711:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
    712:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
    713:   DMDAVecGetArray(properties_da,local_properties,&props);

    715:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
    716:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
    717:   for (ej = sey; ej < sey+my; ej++) {
    718:     for (ei = sex; ei < sex+mx; ei++) {
    719:       /* get coords for the element */
    720:       GetElementCoords(_coords,ei,ej,el_coords);

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

    725:       /* initialise element stiffness matrix */
    726:       PetscMemzero(Ae,sizeof(Ae));
    727:       PetscMemzero(Ge,sizeof(Ge));
    728:       PetscMemzero(De,sizeof(De));
    729:       PetscMemzero(Ce,sizeof(Ce));

    731:       /* form element stiffness matrix */
    732:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
    733:       FormGradientOperatorQ1(Ge,el_coords);
    734:       FormDivergenceOperatorQ1(De,el_coords);
    735:       FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);

    737:       /* insert element matrix into global matrix */
    738:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
    739:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
    740:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
    741:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
    742:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
    743:     }
    744:   }
    745:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    746:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    748:   DMDAVecRestoreArray(cda,coords,&_coords);

    750:   DMDAVecRestoreArray(properties_da,local_properties,&props);
    751:   VecDestroy(&local_properties);
    752:   return(0);
    753: }

    755: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
    756: {
    757:   DM                     cda;
    758:   Vec                    coords;
    759:   DMDACoor2d             **_coords;
    760:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    761:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
    762:   PetscInt               sex,sey,mx,my;
    763:   PetscInt               ei,ej;
    764:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
    765:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
    766:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
    767:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
    768:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    769:   Vec                    local_properties;
    770:   GaussPointCoefficients **props;
    771:   PetscScalar            *prop_eta;
    772:   PetscErrorCode         ierr;

    775:   /* setup for coords */
    776:   DMGetCoordinateDM(stokes_da,&cda);
    777:   DMGetCoordinatesLocal(stokes_da,&coords);
    778:   DMDAVecGetArray(cda,coords,&_coords);

    780:   /* setup for coefficients */
    781:   DMCreateLocalVector(properties_da,&local_properties);
    782:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
    783:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
    784:   DMDAVecGetArray(properties_da,local_properties,&props);

    786:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
    787:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
    788:   for (ej = sey; ej < sey+my; ej++) {
    789:     for (ei = sex; ei < sex+mx; ei++) {
    790:       /* get coords for the element */
    791:       GetElementCoords(_coords,ei,ej,el_coords);

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

    796:       /* initialise element stiffness matrix */
    797:       PetscMemzero(Ae,sizeof(Ae));
    798:       PetscMemzero(Ge,sizeof(Ge));
    799:       PetscMemzero(De,sizeof(De));
    800:       PetscMemzero(Ce,sizeof(Ce));

    802:       /* form element stiffness matrix */
    803:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
    804:       FormGradientOperatorQ1(Ge,el_coords);
    805:       FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);

    807:       /* insert element matrix into global matrix */
    808:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
    809:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
    810:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
    811:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
    812:     }
    813:   }
    814:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    815:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    817:   DMDAVecRestoreArray(cda,coords,&_coords);

    819:   DMDAVecRestoreArray(properties_da,local_properties,&props);
    820:   VecDestroy(&local_properties);
    821:   return(0);
    822: }

    824: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
    825: {
    826:   PetscInt n;

    829:   for (n = 0; n < 4; n++) {
    830:     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];
    831:     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];
    832:     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];
    833:   }
    834:   return(0);
    835: }

    837: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
    838: {
    839:   DM                     cda;
    840:   Vec                    coords;
    841:   DMDACoor2d             **_coords;
    842:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    843:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
    844:   PetscInt               sex,sey,mx,my;
    845:   PetscInt               ei,ej;
    846:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
    847:   PetscScalar            He[NODES_PER_EL*P_DOFS];
    848:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    849:   Vec                    local_properties;
    850:   GaussPointCoefficients **props;
    851:   PetscScalar            *prop_fx,*prop_fy;
    852:   Vec                    local_F;
    853:   StokesDOF              **ff;
    854:   PetscErrorCode         ierr;

    857:   /* setup for coords */
    858:   DMGetCoordinateDM(stokes_da,&cda);
    859:   DMGetCoordinatesLocal(stokes_da,&coords);
    860:   DMDAVecGetArray(cda,coords,&_coords);

    862:   /* setup for coefficients */
    863:   DMGetLocalVector(properties_da,&local_properties);
    864:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
    865:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
    866:   DMDAVecGetArray(properties_da,local_properties,&props);

    868:   /* get acces to the vector */
    869:   DMGetLocalVector(stokes_da,&local_F);
    870:   VecZeroEntries(local_F);
    871:   DMDAVecGetArray(stokes_da,local_F,&ff);

    873:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
    874:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
    875:   for (ej = sey; ej < sey+my; ej++) {
    876:     for (ei = sex; ei < sex+mx; ei++) {
    877:       /* get coords for the element */
    878:       GetElementCoords(_coords,ei,ej,el_coords);

    880:       /* get coefficients for the element */
    881:       prop_fx = props[ej][ei].fx;
    882:       prop_fy = props[ej][ei].fy;

    884:       /* initialise element stiffness matrix */
    885:       PetscMemzero(Fe,sizeof(Fe));
    886:       PetscMemzero(He,sizeof(He));

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

    891:       /* insert element matrix into global matrix */
    892:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);

    894:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
    895:     }
    896:   }

    898:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
    899:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
    900:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
    901:   DMRestoreLocalVector(stokes_da,&local_F);

    903:   DMDAVecRestoreArray(cda,coords,&_coords);

    905:   DMDAVecRestoreArray(properties_da,local_properties,&props);
    906:   DMRestoreLocalVector(properties_da,&local_properties);
    907:   return(0);
    908: }

    910: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
    911: {
    912:   DM             da,cda;
    913:   Vec            X;
    914:   StokesDOF      **_stokes;
    915:   Vec            coords;
    916:   DMDACoor2d     **_coords;
    917:   PetscInt       si,sj,ei,ej,i,j;

    921:   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);
    922:   DMSetFromOptions(da);
    923:   DMSetUp(da);
    924:   DMDASetFieldName(da,0,"anlytic_Vx");
    925:   DMDASetFieldName(da,1,"anlytic_Vy");
    926:   DMDASetFieldName(da,2,"analytic_P");

    928:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);

    930:   DMGetCoordinatesLocal(da,&coords);
    931:   DMGetCoordinateDM(da,&cda);
    932:   DMDAVecGetArray(cda,coords,&_coords);

    934:   DMCreateGlobalVector(da,&X);
    935:   DMDAVecGetArray(da,X,&_stokes);

    937:   DMDAGetCorners(da,&si,&sj,0,&ei,&ej,0);
    938:   for (j = sj; j < sj+ej; j++) {
    939:     for (i = si; i < si+ei; i++) {
    940:       PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];

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

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

    947:       _stokes[j][i].u_dof = vel[0];
    948:       _stokes[j][i].v_dof = vel[1];
    949:       _stokes[j][i].p_dof = pressure;
    950:     }
    951:   }
    952:   DMDAVecRestoreArray(da,X,&_stokes);
    953:   DMDAVecRestoreArray(cda,coords,&_coords);

    955:   *_da = da;
    956:   *_X  = X;
    957:   return(0);
    958: }

    960: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
    961: {
    963:   /* get the nodal fields */
    964:   nodal_fields[0].u_dof = fields[ej][ei].u_dof;     nodal_fields[0].v_dof = fields[ej][ei].v_dof;     nodal_fields[0].p_dof = fields[ej][ei].p_dof;
    965:   nodal_fields[1].u_dof = fields[ej+1][ei].u_dof;   nodal_fields[1].v_dof = fields[ej+1][ei].v_dof;   nodal_fields[1].p_dof = fields[ej+1][ei].p_dof;
    966:   nodal_fields[2].u_dof = fields[ej+1][ei+1].u_dof; nodal_fields[2].v_dof = fields[ej+1][ei+1].v_dof; nodal_fields[2].p_dof = fields[ej+1][ei+1].p_dof;
    967:   nodal_fields[3].u_dof = fields[ej][ei+1].u_dof;   nodal_fields[3].v_dof = fields[ej][ei+1].v_dof;   nodal_fields[3].p_dof = fields[ej][ei+1].p_dof;
    968:   return(0);
    969: }

    971: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
    972: {
    973:   DM             cda;
    974:   Vec            coords,X_analytic_local,X_local;
    975:   DMDACoor2d     **_coords;
    976:   PetscInt       sex,sey,mx,my;
    977:   PetscInt       ei,ej;
    978:   PetscScalar    el_coords[NODES_PER_EL*NSD];
    979:   StokesDOF      **stokes_analytic,**stokes;
    980:   StokesDOF      stokes_analytic_e[4],stokes_e[4];

    982:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    983:   PetscScalar    Ni_p[NODES_PER_EL];
    984:   PetscInt       ngp;
    985:   PetscScalar    gp_xi[GAUSS_POINTS][2];
    986:   PetscScalar    gp_weight[GAUSS_POINTS];
    987:   PetscInt       p,i;
    988:   PetscScalar    J_p,fac;
    989:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
    990:   PetscInt       M;
    991:   PetscReal      xymin[2],xymax[2];

    995:   /* define quadrature rule */
    996:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    998:   /* setup for coords */
    999:   DMGetCoordinateDM(stokes_da,&cda);
    1000:   DMGetCoordinatesLocal(stokes_da,&coords);
    1001:   DMDAVecGetArray(cda,coords,&_coords);

    1003:   /* setup for analytic */
    1004:   DMCreateLocalVector(stokes_da,&X_analytic_local);
    1005:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
    1006:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
    1007:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

    1009:   /* setup for solution */
    1010:   DMCreateLocalVector(stokes_da,&X_local);
    1011:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
    1012:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
    1013:   DMDAVecGetArray(stokes_da,X_local,&stokes);

    1015:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
    1016:   DMGetBoundingBox(stokes_da,xymin,xymax);

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

    1020:   tp_L2 = tu_L2 = tu_H1 = 0.0;

    1022:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
    1023:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
    1024:   for (ej = sey; ej < sey+my; ej++) {
    1025:     for (ei = sex; ei < sex+mx; ei++) {
    1026:       /* get coords for the element */
    1027:       GetElementCoords(_coords,ei,ej,el_coords);
    1028:       StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
    1029:       StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);

    1031:       /* evaluate integral */
    1032:       p_e_L2 = 0.0;
    1033:       u_e_L2 = 0.0;
    1034:       u_e_H1 = 0.0;
    1035:       for (p = 0; p < ngp; p++) {
    1036:         ConstructQ12D_Ni(gp_xi[p],Ni_p);
    1037:         ConstructQ12D_GNi(gp_xi[p],GNi_p);
    1038:         ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
    1039:         fac = gp_weight[p]*J_p;

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

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

    1046:           u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
    1047:           v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
    1048:           u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);

    1050:           u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
    1051:                                +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error               /* du/dy */
    1052:                                +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
    1053:                                +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error);             /* dv/dy */
    1054:         }
    1055:       }

    1057:       tp_L2 += p_e_L2;
    1058:       tu_L2 += u_e_L2;
    1059:       tu_H1 += u_e_H1;
    1060:     }
    1061:   }
    1062:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
    1063:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
    1064:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
    1065:   p_L2 = PetscSqrtScalar(p_L2);
    1066:   u_L2 = PetscSqrtScalar(u_L2);
    1067:   u_H1 = PetscSqrtScalar(u_H1);

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

    1071:   DMDAVecRestoreArray(cda,coords,&_coords);

    1073:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
    1074:   VecDestroy(&X_analytic_local);
    1075:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
    1076:   VecDestroy(&X_local);
    1077:   return(0);
    1078: }

    1080: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
    1081: {
    1082:   DM                     da_Stokes,da_prop;
    1083:   PetscInt               u_dof,p_dof,dof,stencil_width;
    1084:   Mat                    A,B;
    1085:   DM                     prop_cda,vel_cda;
    1086:   Vec                    prop_coords,vel_coords;
    1087:   PetscInt               si,sj,nx,ny,i,j,p;
    1088:   Vec                    f,X;
    1089:   PetscInt               prop_dof,prop_stencil_width;
    1090:   Vec                    properties,l_properties;
    1091:   PetscReal              dx,dy;
    1092:   PetscInt               M,N;
    1093:   DMDACoor2d             **_prop_coords,**_vel_coords;
    1094:   GaussPointCoefficients **element_props;
    1095:   PetscInt               its;
    1096:   KSP                    ksp_S;
    1097:   PetscInt               coefficient_structure = 0;
    1098:   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
    1099:   PetscBool              use_gp_coords = PETSC_FALSE,set,output_gnuplot = PETSC_FALSE,glvis = PETSC_FALSE;
    1100:   char                   filename[PETSC_MAX_PATH_LEN];
    1101:   PetscErrorCode         ierr;

    1105:   PetscOptionsGetBool(NULL,NULL,"-gnuplot",&output_gnuplot,NULL);
    1106:   PetscOptionsGetBool(NULL,NULL,"-glvis",&glvis,NULL);

    1108:   /* Generate the da for velocity and pressure */
    1109:   /*
    1110:   We use Q1 elements for the temperature.
    1111:   FEM has a 9-point stencil (BOX) or connectivity pattern
    1112:   Num nodes in each direction is mx+1, my+1
    1113:   */
    1114:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
    1115:   p_dof         = P_DOFS; /* p - pressure */
    1116:   dof           = u_dof+p_dof;
    1117:   stencil_width = 1;
    1118:   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);

    1120:   DMSetMatType(da_Stokes,MATAIJ);
    1121:   DMSetFromOptions(da_Stokes);
    1122:   DMSetUp(da_Stokes);
    1123:   DMDASetFieldName(da_Stokes,0,"Vx");
    1124:   DMDASetFieldName(da_Stokes,1,"Vy");
    1125:   DMDASetFieldName(da_Stokes,2,"P");

    1127:   /* unit box [0,1] x [0,1] */
    1128:   DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);

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

    1135:   prop_dof           = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
    1136:   prop_stencil_width = 0;
    1137:   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);
    1138:   DMSetFromOptions(da_prop);
    1139:   DMSetUp(da_prop);
    1140:   PetscFree(lx);
    1141:   PetscFree(ly);

    1143:   /* define centroid positions */
    1144:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
    1145:   dx   = 1.0/((PetscReal)(M));
    1146:   dy   = 1.0/((PetscReal)(N));

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

    1150:   /* define coefficients */
    1151:   PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);

    1153:   DMCreateGlobalVector(da_prop,&properties);
    1154:   DMCreateLocalVector(da_prop,&l_properties);
    1155:   DMDAVecGetArray(da_prop,l_properties,&element_props);

    1157:   DMGetCoordinateDM(da_prop,&prop_cda);
    1158:   DMGetCoordinatesLocal(da_prop,&prop_coords);
    1159:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

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

    1163:   DMGetCoordinateDM(da_Stokes,&vel_cda);
    1164:   DMGetCoordinatesLocal(da_Stokes,&vel_coords);
    1165:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);

    1167:   /* interpolate the coordinates */
    1168:   for (j = sj; j < sj+ny; j++) {
    1169:     for (i = si; i < si+nx; i++) {
    1170:       PetscInt    ngp;
    1171:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
    1172:       PetscScalar el_coords[8];

    1174:       GetElementCoords(_vel_coords,i,j,el_coords);
    1175:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    1177:       for (p = 0; p < GAUSS_POINTS; p++) {
    1178:         PetscScalar gp_x,gp_y;
    1179:         PetscInt    n;
    1180:         PetscScalar xi_p[2],Ni_p[4];

    1182:         xi_p[0] = gp_xi[p][0];
    1183:         xi_p[1] = gp_xi[p][1];
    1184:         ConstructQ12D_Ni(xi_p,Ni_p);

    1186:         gp_x = 0.0;
    1187:         gp_y = 0.0;
    1188:         for (n = 0; n < NODES_PER_EL; n++) {
    1189:           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
    1190:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
    1191:         }
    1192:         element_props[j][i].gp_coords[2*p]   = gp_x;
    1193:         element_props[j][i].gp_coords[2*p+1] = gp_y;
    1194:       }
    1195:     }
    1196:   }

    1198:   /* define the coefficients */
    1199:   PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,NULL);

    1201:   for (j = sj; j < sj+ny; j++) {
    1202:     for (i = si; i < si+nx; i++) {
    1203:       PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
    1204:       PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
    1205:       PetscReal coord_x,coord_y;

    1207:       if (coefficient_structure == 0) {
    1208:         PetscReal opts_eta0,opts_eta1,opts_xc;
    1209:         PetscInt  opts_nz;

    1211:         opts_eta0 = 1.0;
    1212:         opts_eta1 = 1.0;
    1213:         opts_xc   = 0.5;
    1214:         opts_nz   = 1;

    1216:         PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
    1217:         PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
    1218:         PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
    1219:         PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);

    1221:         for (p = 0; p < GAUSS_POINTS; p++) {
    1222:           coord_x = centroid_x;
    1223:           coord_y = centroid_y;
    1224:           if (use_gp_coords) {
    1225:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
    1226:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
    1227:           }

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

    1232:           element_props[j][i].fx[p] = 0.0;
    1233:           element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
    1234:         }
    1235:       } else if (coefficient_structure == 1) { /* square sinker */
    1236:         PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;

    1238:         opts_eta0 = 1.0;
    1239:         opts_eta1 = 1.0;
    1240:         opts_dx   = 0.50;
    1241:         opts_dy   = 0.50;

    1243:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
    1244:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
    1245:         PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
    1246:         PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);

    1248:         for (p = 0; p < GAUSS_POINTS; p++) {
    1249:           coord_x = centroid_x;
    1250:           coord_y = centroid_y;
    1251:           if (use_gp_coords) {
    1252:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
    1253:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
    1254:           }

    1256:           element_props[j][i].eta[p] = opts_eta0;
    1257:           element_props[j][i].fx[p]  = 0.0;
    1258:           element_props[j][i].fy[p]  = 0.0;

    1260:           if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
    1261:             if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
    1262:               element_props[j][i].eta[p] =  opts_eta1;
    1263:               element_props[j][i].fx[p]  =  0.0;
    1264:               element_props[j][i].fy[p]  = -1.0;
    1265:             }
    1266:           }
    1267:         }
    1268:       } else if (coefficient_structure == 2) { /* circular sinker */
    1269:         PetscReal opts_eta0,opts_eta1,opts_r,radius2;

    1271:         opts_eta0 = 1.0;
    1272:         opts_eta1 = 1.0;
    1273:         opts_r    = 0.25;

    1275:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
    1276:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
    1277:         PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);

    1279:         for (p = 0; p < GAUSS_POINTS; p++) {
    1280:           coord_x = centroid_x;
    1281:           coord_y = centroid_y;
    1282:           if (use_gp_coords) {
    1283:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
    1284:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
    1285:           }

    1287:           element_props[j][i].eta[p] = opts_eta0;
    1288:           element_props[j][i].fx[p]  = 0.0;
    1289:           element_props[j][i].fy[p]  = 0.0;

    1291:           radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
    1292:           if (radius2 < opts_r*opts_r) {
    1293:             element_props[j][i].eta[p] =  opts_eta1;
    1294:             element_props[j][i].fx[p]  =  0.0;
    1295:             element_props[j][i].fy[p]  = -1.0;
    1296:           }
    1297:         }
    1298:       } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
    1299:         PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;

    1301:         opts_eta0 = 1.0;
    1302:         opts_eta1 = 1.0;
    1303:         opts_r    = 0.25;
    1304:         opts_c0x  = 0.35;       /* circle center */
    1305:         opts_c0y  = 0.35;
    1306:         opts_s0x  = 0.7;       /* square center */
    1307:         opts_s0y  = 0.7;
    1308:         opts_dx   = 0.25;
    1309:         opts_dy   = 0.25;
    1310:         opts_phi  = 25;

    1312:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
    1313:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
    1314:         PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
    1315:         PetscOptionsGetReal(NULL,NULL,"-sinker_c0x",&opts_c0x,NULL);
    1316:         PetscOptionsGetReal(NULL,NULL,"-sinker_c0y",&opts_c0y,NULL);
    1317:         PetscOptionsGetReal(NULL,NULL,"-sinker_s0x",&opts_s0x,NULL);
    1318:         PetscOptionsGetReal(NULL,NULL,"-sinker_s0y",&opts_s0y,NULL);
    1319:         PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
    1320:         PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
    1321:         PetscOptionsGetReal(NULL,NULL,"-sinker_phi",&opts_phi,NULL);
    1322:         opts_phi *= PETSC_PI / 180;

    1324:         for (p = 0; p < GAUSS_POINTS; p++) {
    1325:           coord_x = centroid_x;
    1326:           coord_y = centroid_y;
    1327:           if (use_gp_coords) {
    1328:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
    1329:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
    1330:           }

    1332:           element_props[j][i].eta[p] = opts_eta0;
    1333:           element_props[j][i].fx[p]  = 0.0;
    1334:           element_props[j][i].fy[p]  = -0.2;

    1336:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
    1337:           if (radius2 < opts_r*opts_r
    1338:               || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
    1339:                   && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
    1340:             element_props[j][i].eta[p] =  opts_eta1;
    1341:             element_props[j][i].fx[p]  =  0.0;
    1342:             element_props[j][i].fy[p]  = -1.0;
    1343:           }
    1344:         }
    1345:       } else if (coefficient_structure == 4) { /* subdomain jump */
    1346:         PetscReal opts_mag,opts_eta0;
    1347:         PetscInt  opts_nz,px,py;
    1348:         PetscBool jump;

    1350:         opts_mag  = 1.0;
    1351:         opts_eta0 = 1.0;
    1352:         opts_nz   = 1;

    1354:         PetscOptionsGetReal(NULL,NULL,"-jump_eta0",&opts_eta0,NULL);
    1355:         PetscOptionsGetReal(NULL,NULL,"-jump_magnitude",&opts_mag,NULL);
    1356:         PetscOptionsGetInt(NULL,NULL,"-jump_nz",&opts_nz,NULL);
    1357:         DMDAGetInfo(da_Stokes,NULL,NULL,NULL,NULL,&px,&py,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
    1358:         if (px%2) {
    1359:           jump = (PetscBool)(PetscGlobalRank%2);
    1360:         } else {
    1361:           jump = (PetscBool)((PetscGlobalRank/px)%2 ? PetscGlobalRank%2 : !(PetscGlobalRank%2));
    1362:         }
    1363:         for (p = 0; p < GAUSS_POINTS; p++) {
    1364:           coord_x = centroid_x;
    1365:           coord_y = centroid_y;
    1366:           if (use_gp_coords) {
    1367:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
    1368:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
    1369:           }

    1371:           element_props[j][i].eta[p] = jump ? PetscPowReal(10.0,opts_mag) : opts_eta0;
    1372:           element_props[j][i].fx[p]  = 0.0;
    1373:           element_props[j][i].fy[p]  = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
    1374:         }
    1375:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
    1376:     }
    1377:   }
    1378:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

    1380:   DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);

    1382:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
    1383:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
    1384:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);

    1386:   if (output_gnuplot) {
    1387:     DMDACoordViewGnuplot2d(da_Stokes,"mesh");
    1388:     DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coefficients for Stokes eqn.","properties");
    1389:   }

    1391:   if (glvis) {
    1392:     Vec         glv_prop,etaf;
    1393:     PetscViewer view;
    1394:     PetscInt    dim;
    1395:     const char  *fec = {"FiniteElementCollection: L2_2D_P1"};

    1397:     DMGetDimension(da_Stokes,&dim);
    1398:     VecCreateSeq(PETSC_COMM_SELF,GAUSS_POINTS*mx*mx,&etaf);
    1399:     PetscObjectSetName((PetscObject)etaf,"viscosity");
    1400:     PetscViewerGLVisOpen(PETSC_COMM_WORLD,PETSC_VIEWER_GLVIS_SOCKET,NULL,PETSC_DECIDE,&view);
    1401:     PetscViewerGLVisSetFields(view,1,&fec,&dim,glvis_extract_eta,(PetscObject*)&etaf,da_prop,NULL);
    1402:     DMGetLocalVector(da_prop,&glv_prop);
    1403:     DMGlobalToLocalBegin(da_prop,properties,INSERT_VALUES,glv_prop);
    1404:     DMGlobalToLocalEnd(da_prop,properties,INSERT_VALUES,glv_prop);
    1405:     VecSetDM(etaf,da_Stokes);
    1406:     VecView(glv_prop,view);
    1407:     PetscViewerDestroy(&view);
    1408:     DMRestoreLocalVector(da_prop,&glv_prop);
    1409:     VecDestroy(&etaf);
    1410:   }

    1412:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
    1413:   DMCreateMatrix(da_Stokes,&A);
    1414:   DMCreateMatrix(da_Stokes,&B);
    1415:   DMCreateGlobalVector(da_Stokes,&f);
    1416:   DMCreateGlobalVector(da_Stokes,&X);

    1418:   /* assemble A11 */
    1419:   MatZeroEntries(A);
    1420:   MatZeroEntries(B);
    1421:   VecZeroEntries(f);

    1423:   AssembleA_Stokes(A,da_Stokes,da_prop,properties);
    1424:   AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
    1425:   /* build force vector */
    1426:   AssembleF_Stokes(f,da_Stokes,da_prop,properties);

    1428:   DMDABCApplyFreeSlip(da_Stokes,A,f);
    1429:   DMDABCApplyFreeSlip(da_Stokes,B,NULL);

    1431:   /* SOLVE */
    1432:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
    1433:   KSPSetOptionsPrefix(ksp_S,"stokes_");
    1434:   KSPSetDM(ksp_S,da_Stokes);
    1435:   KSPSetDMActive(ksp_S,PETSC_FALSE);
    1436:   KSPSetOperators(ksp_S,A,B);
    1437:   KSPSetFromOptions(ksp_S);
    1438:   {
    1439:     PC             pc;
    1440:     const PetscInt ufields[] = {0,1},pfields[1] = {2};

    1442:     KSPGetPC(ksp_S,&pc);
    1443:     PCFieldSplitSetBlockSize(pc,3);
    1444:     PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
    1445:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
    1446:   }

    1448:   {
    1449:     PC        pc;
    1450:     PetscBool same = PETSC_FALSE;
    1451:     KSPGetPC(ksp_S,&pc);
    1452:     PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&same);
    1453:     if (same) {
    1454:       PetscBool usedivmat = PETSC_FALSE;
    1455:       KSPSetOperators(ksp_S,A,A);

    1457:       PetscOptionsGetBool(NULL,NULL,"-stokes_pc_bddc_use_divergence",&usedivmat,NULL);
    1458:       if (usedivmat) {
    1459:         IS       *fields,vel;
    1460:         PetscInt i,nf;

    1462:         DMCreateFieldDecomposition(da_Stokes,&nf,NULL,&fields,NULL);
    1463:         ISConcatenate(PETSC_COMM_WORLD,2,fields,&vel);
    1464:         MatZeroRowsIS(B,fields[2],1.0,NULL,NULL); /* we put 1.0 on the diagonal to pick the pressure average too */
    1465:         MatTranspose(B,MAT_INPLACE_MATRIX,&B);
    1466:         MatZeroRowsIS(B,vel,0.0,NULL,NULL);
    1467:         ISDestroy(&vel);
    1468:         PCBDDCSetDivergenceMat(pc,B,PETSC_FALSE,NULL);
    1469:         for (i=0;i<nf;i++) {
    1470:           ISDestroy(&fields[i]);
    1471:         }
    1472:         PetscFree(fields);
    1473:       }
    1474:     }
    1475:   }

    1477:   {
    1478:     PC        pc_S;
    1479:     KSP       *sub_ksp,ksp_U;
    1480:     PetscInt  nsplits;
    1481:     DM        da_U;
    1482:     PetscBool is_pcfs;

    1484:     KSPSetUp(ksp_S);
    1485:     KSPGetPC(ksp_S,&pc_S);

    1487:     is_pcfs = PETSC_FALSE;
    1488:     PetscObjectTypeCompare((PetscObject)pc_S,PCFIELDSPLIT,&is_pcfs);

    1490:     if (is_pcfs) {
    1491:       PCFieldSplitGetSubKSP(pc_S,&nsplits,&sub_ksp);
    1492:       ksp_U = sub_ksp[0];
    1493:       PetscFree(sub_ksp);

    1495:       if (nsplits == 2) {
    1496:         DMDACreateCompatibleDMDA(da_Stokes,2,&da_U);

    1498:         KSPSetDM(ksp_U,da_U);
    1499:         KSPSetDMActive(ksp_U,PETSC_FALSE);
    1500:         DMDestroy(&da_U);
    1501:       }
    1502:     }
    1503:   }

    1505:   KSPSolve(ksp_S,f,X);

    1507:   PetscOptionsGetString(NULL,NULL,"-o",filename,sizeof(filename),&set);
    1508:   if (set) {
    1509:     char        *ext;
    1510:     PetscViewer viewer;
    1511:     PetscBool   flg;
    1512:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
    1513:     PetscStrrchr(filename,'.',&ext);
    1514:     PetscStrcmp("vts",ext,&flg);
    1515:     if (flg) {
    1516:       PetscViewerSetType(viewer,PETSCVIEWERVTK);
    1517:     } else {
    1518:       PetscViewerSetType(viewer,PETSCVIEWERBINARY);
    1519:     }
    1520:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
    1521:     PetscViewerFileSetName(viewer,filename);
    1522:     VecView(X,viewer);
    1523:     PetscViewerDestroy(&viewer);
    1524:   }
    1525:   if (output_gnuplot) {
    1526:     DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
    1527:   }

    1529:   if (glvis) {
    1530:     PetscViewer view;

    1532:     PetscViewerCreate(PETSC_COMM_WORLD,&view);
    1533:     PetscViewerSetType(view,PETSCVIEWERGLVIS);
    1534:     VecView(X,view);
    1535:     PetscViewerDestroy(&view);
    1536:   }

    1538:   KSPGetIterationNumber(ksp_S,&its);

    1540:   if (coefficient_structure == 0) {
    1541:     PetscReal opts_eta0,opts_eta1,opts_xc;
    1542:     PetscInt  opts_nz,N;
    1543:     DM        da_Stokes_analytic;
    1544:     Vec       X_analytic;
    1545:     PetscReal nrm1[3],nrm2[3],nrmI[3];

    1547:     opts_eta0 = 1.0;
    1548:     opts_eta1 = 1.0;
    1549:     opts_xc   = 0.5;
    1550:     opts_nz   = 1;

    1552:     PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
    1553:     PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
    1554:     PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
    1555:     PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);

    1557:     DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
    1558:     if (output_gnuplot) {
    1559:       DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
    1560:     }
    1561:     DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);

    1563:     VecAXPY(X_analytic,-1.0,X);
    1564:     VecGetSize(X_analytic,&N);
    1565:     N    = N/3;

    1567:     VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
    1568:     VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
    1569:     VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);

    1571:     VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
    1572:     VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
    1573:     VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);

    1575:     VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
    1576:     VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
    1577:     VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);

    1579:     DMDestroy(&da_Stokes_analytic);
    1580:     VecDestroy(&X_analytic);
    1581:   }

    1583:   KSPDestroy(&ksp_S);
    1584:   VecDestroy(&X);
    1585:   VecDestroy(&f);
    1586:   MatDestroy(&A);
    1587:   MatDestroy(&B);

    1589:   DMDestroy(&da_Stokes);
    1590:   DMDestroy(&da_prop);

    1592:   VecDestroy(&properties);
    1593:   VecDestroy(&l_properties);
    1594:   return(0);
    1595: }

    1597: int main(int argc,char **args)
    1598: {
    1600:   PetscInt       mx,my;

    1602:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
    1603:   mx   = my = 10;
    1604:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
    1605:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
    1606:   solve_stokes_2d_coupled(mx,my);
    1607:   PetscFinalize();
    1608:   return ierr;
    1609: }

    1611: /* -------------------------- helpers for boundary conditions -------------------------------- */
    1612: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
    1613: {
    1614:   DM                     cda;
    1615:   Vec                    coords;
    1616:   PetscInt               si,sj,nx,ny,i,j;
    1617:   PetscInt               M,N;
    1618:   DMDACoor2d             **_coords;
    1619:   const PetscInt         *g_idx;
    1620:   PetscInt               *bc_global_ids;
    1621:   PetscScalar            *bc_vals;
    1622:   PetscInt               nbcs;
    1623:   PetscInt               n_dofs;
    1624:   PetscErrorCode         ierr;
    1625:   ISLocalToGlobalMapping ltogm;

    1628:   DMGetLocalToGlobalMapping(da,&ltogm);
    1629:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1631:   DMGetCoordinateDM(da,&cda);
    1632:   DMGetCoordinatesLocal(da,&coords);
    1633:   DMDAVecGetArray(cda,coords,&_coords);
    1634:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1635:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1637:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
    1638:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

    1643:   i = nx-1;
    1644:   for (j = 0; j < ny; j++) {
    1645:     PetscInt local_id;

    1647:     local_id = i+j*nx;

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

    1651:     bc_vals[j] =  0.0;
    1652:   }
    1653:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1654:   nbcs = 0;
    1655:   if ((si+nx) == (M)) nbcs = ny;

    1657:   if (b) {
    1658:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1659:     VecAssemblyBegin(b);
    1660:     VecAssemblyEnd(b);
    1661:   }
    1662:   if (A) {
    1663:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1664:   }

    1666:   PetscFree(bc_vals);
    1667:   PetscFree(bc_global_ids);

    1669:   DMDAVecRestoreArray(cda,coords,&_coords);
    1670:   return(0);
    1671: }

    1673: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
    1674: {
    1675:   DM                     cda;
    1676:   Vec                    coords;
    1677:   PetscInt               si,sj,nx,ny,i,j;
    1678:   PetscInt               M,N;
    1679:   DMDACoor2d             **_coords;
    1680:   const PetscInt         *g_idx;
    1681:   PetscInt               *bc_global_ids;
    1682:   PetscScalar            *bc_vals;
    1683:   PetscInt               nbcs;
    1684:   PetscInt               n_dofs;
    1685:   PetscErrorCode         ierr;
    1686:   ISLocalToGlobalMapping ltogm;

    1689:   DMGetLocalToGlobalMapping(da,&ltogm);
    1690:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1692:   DMGetCoordinateDM(da,&cda);
    1693:   DMGetCoordinatesLocal(da,&coords);
    1694:   DMDAVecGetArray(cda,coords,&_coords);
    1695:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1696:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1698:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
    1699:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

    1704:   i = 0;
    1705:   for (j = 0; j < ny; j++) {
    1706:     PetscInt local_id;

    1708:     local_id = i+j*nx;

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

    1712:     bc_vals[j] =  0.0;
    1713:   }
    1714:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1715:   nbcs = 0;
    1716:   if (si == 0) nbcs = ny;

    1718:   if (b) {
    1719:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1720:     VecAssemblyBegin(b);
    1721:     VecAssemblyEnd(b);
    1722:   }

    1724:   if (A) {
    1725:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1726:   }

    1728:   PetscFree(bc_vals);
    1729:   PetscFree(bc_global_ids);

    1731:   DMDAVecRestoreArray(cda,coords,&_coords);
    1732:   return(0);
    1733: }

    1735: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
    1736: {
    1737:   DM                     cda;
    1738:   Vec                    coords;
    1739:   PetscInt               si,sj,nx,ny,i,j;
    1740:   PetscInt               M,N;
    1741:   DMDACoor2d             **_coords;
    1742:   const PetscInt         *g_idx;
    1743:   PetscInt               *bc_global_ids;
    1744:   PetscScalar            *bc_vals;
    1745:   PetscInt               nbcs;
    1746:   PetscInt               n_dofs;
    1747:   PetscErrorCode         ierr;
    1748:   ISLocalToGlobalMapping ltogm;

    1751:   DMGetLocalToGlobalMapping(da,&ltogm);
    1752:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1754:   DMGetCoordinateDM(da,&cda);
    1755:   DMGetCoordinatesLocal(da,&coords);
    1756:   DMDAVecGetArray(cda,coords,&_coords);
    1757:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1758:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1760:   PetscMalloc1(nx,&bc_global_ids);
    1761:   PetscMalloc1(nx,&bc_vals);

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

    1766:   j = ny-1;
    1767:   for (i = 0; i < nx; i++) {
    1768:     PetscInt local_id;

    1770:     local_id = i+j*nx;

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

    1774:     bc_vals[i] =  0.0;
    1775:   }
    1776:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1777:   nbcs = 0;
    1778:   if ((sj+ny) == (N)) nbcs = nx;

    1780:   if (b) {
    1781:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1782:     VecAssemblyBegin(b);
    1783:     VecAssemblyEnd(b);
    1784:   }
    1785:   if (A) {
    1786:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1787:   }

    1789:   PetscFree(bc_vals);
    1790:   PetscFree(bc_global_ids);

    1792:   DMDAVecRestoreArray(cda,coords,&_coords);
    1793:   return(0);
    1794: }

    1796: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
    1797: {
    1798:   DM                     cda;
    1799:   Vec                    coords;
    1800:   PetscInt               si,sj,nx,ny,i,j;
    1801:   PetscInt               M,N;
    1802:   DMDACoor2d             **_coords;
    1803:   const PetscInt         *g_idx;
    1804:   PetscInt               *bc_global_ids;
    1805:   PetscScalar            *bc_vals;
    1806:   PetscInt               nbcs;
    1807:   PetscInt               n_dofs;
    1808:   PetscErrorCode         ierr;
    1809:   ISLocalToGlobalMapping ltogm;

    1812:   DMGetLocalToGlobalMapping(da,&ltogm);
    1813:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1815:   DMGetCoordinateDM(da,&cda);
    1816:   DMGetCoordinatesLocal(da,&coords);
    1817:   DMDAVecGetArray(cda,coords,&_coords);
    1818:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1819:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1821:   PetscMalloc1(nx,&bc_global_ids);
    1822:   PetscMalloc1(nx,&bc_vals);

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

    1827:   j = 0;
    1828:   for (i = 0; i < nx; i++) {
    1829:     PetscInt local_id;

    1831:     local_id = i+j*nx;

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

    1835:     bc_vals[i] =  0.0;
    1836:   }
    1837:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1838:   nbcs = 0;
    1839:   if (sj == 0) nbcs = nx;

    1841:   if (b) {
    1842:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1843:     VecAssemblyBegin(b);
    1844:     VecAssemblyEnd(b);
    1845:   }
    1846:   if (A) {
    1847:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1848:   }

    1850:   PetscFree(bc_vals);
    1851:   PetscFree(bc_global_ids);

    1853:   DMDAVecRestoreArray(cda,coords,&_coords);
    1854:   return(0);
    1855: }

    1857: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
    1858: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
    1859: {

    1863:   BCApplyZero_NORTH(da_Stokes,1,A,f);
    1864:   BCApplyZero_EAST(da_Stokes,0,A,f);
    1865:   BCApplyZero_SOUTH(da_Stokes,1,A,f);
    1866:   BCApplyZero_WEST(da_Stokes,0,A,f);
    1867:   return(0);
    1868: }

    1871: /*TEST

    1873:    build:
    1874:       requires: !complex !single

    1876:    test:
    1877:       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

    1879:    test:
    1880:       suffix: 2
    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_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
    1882:       output_file: output/ex43_1.out

    1884:    test:
    1885:       suffix: 3
    1886:       nsize: 4
    1887:       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

    1889:    test:
    1890:       suffix: 4
    1891:       nsize: 4
    1892:       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

    1894:    test:
    1895:       suffix: 5
    1896:       nsize: 4
    1897:       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

    1899:    test:
    1900:       suffix: 6
    1901:       nsize: 8
    1902:       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

    1904:    test:
    1905:       suffix: bjacobi
    1906:       nsize: 4
    1907:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason

    1909:    test:
    1910:       suffix: bjacobi_baij
    1911:       nsize: 4
    1912:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
    1913:       output_file: output/ex43_bjacobi.out

    1915:    test:
    1916:       suffix: nested_gmg
    1917:       nsize: 4
    1918:       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

    1920:    test:
    1921:       suffix: fetidp
    1922:       nsize: 8
    1923:       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

    1925:    test:
    1926:       suffix: fetidp_unsym
    1927:       nsize: 8
    1928:       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

    1930:    test:
    1931:       suffix: bddc_stokes_deluxe
    1932:       nsize: 8
    1933:       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

    1935:    test:
    1936:       suffix: bddc_stokes_subdomainjump_deluxe
    1937:       nsize: 9
    1938:       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

    1941: TEST*/