Actual source code: ex70.c

      1: static char help[] = "------------------------------------------------------------------------------------------------------------------------------ \n\
      2:   Solves the time-dependent incompressible, variable viscosity Stokes equation in 2D driven by buouyancy variations. \n\
      3:   Time-dependence is introduced by evolving the density (rho) and viscosity (eta) according to \n\
      4:     D \\rho / Dt = 0    and    D \\eta / Dt = 0 \n\
      5:   The Stokes problem is discretized using Q1-Q1 finite elements, stabilized with Bochev's polynomial projection method. \n\
      6:   The hyperbolic evolution equation for density is discretized using a variant of the Particle-In-Cell (PIC) method. \n\
      7:   The DMDA object is used to define the FE problem, whilst DMSwarm provides support for the PIC method. \n\
      8:   Material points (particles) store density and viscosity. The particles are advected with the fluid velocity using RK1. \n\
      9:   At each time step, the value of density and viscosity stored on each particle are projected into a Q1 function space \n\
     10:   and then interpolated onto the Gauss quadrature points. \n\
     11:   The model problem defined in this example is the iso-viscous Rayleigh-Taylor instability (case 1a) from: \n\
     12:     \"A comparison of methods for the modeling of thermochemical convection\" \n\
     13:     P.E. van Keken, S.D. King, H. Schmeling, U.R. Christensen, D. Neumeister and M.-P. Doin, \n\
     14:     Journal of Geophysical Research, vol 102 (B10), 477--499 (1997) \n\
     15:   Note that whilst the model problem defined is for an iso-viscous, the implementation in this example supports \n\
     16:   variable viscoity formulations. \n\
     17:   This example is based on src/ksp/ksp/tutorials/ex43.c \n\
     18:   Options: \n\
     19:     -mx        : Number of elements in the x-direction \n\
     20:     -my        : Number of elements in the y-direction \n\
     21:     -mxy       : Number of elements in the x- and y-directions \n\
     22:     -nt        : Number of time steps \n\
     23:     -dump_freq : Frequency of output file creation \n\
     24:     -ppcell    : Number of times the reference cell is sub-divided \n\
     25:     -randomize_coords : Apply a random shift to each particle coordinate in the range [-fac*dh,0.fac*dh] \n\
     26:     -randomize_fac    : Set the scaling factor for the random shift (default = 0.25)\n";

     28: /* Contributed by Dave May */

     30: #include <petscksp.h>
     31: #include <petscdm.h>
     32: #include <petscdmda.h>
     33: #include <petscdmswarm.h>
     34: #include <petsc/private/dmimpl.h>

     36: static PetscErrorCode DMDAApplyBoundaryConditions(DM,Mat,Vec);

     38: #define NSD            2 /* number of spatial dimensions */
     39: #define NODES_PER_EL   4 /* nodes per element */
     40: #define U_DOFS         2 /* degrees of freedom per velocity node */
     41: #define P_DOFS         1 /* degrees of freedom per pressure node */
     42: #define GAUSS_POINTS   4

     45: static void EvaluateBasis_Q1(PetscScalar _xi[],PetscScalar N[])
     46: {
     47:   PetscScalar xi  = _xi[0];
     48:   PetscScalar eta = _xi[1];

     50:   N[0] = 0.25*(1.0-xi)*(1.0-eta);
     51:   N[1] = 0.25*(1.0+xi)*(1.0-eta);
     52:   N[2] = 0.25*(1.0+xi)*(1.0+eta);
     53:   N[3] = 0.25*(1.0-xi)*(1.0+eta);
     54: }

     56: static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[],PetscScalar dN[][NODES_PER_EL])
     57: {
     58:   PetscScalar xi  = _xi[0];
     59:   PetscScalar eta = _xi[1];

     61:   dN[0][0] = -0.25*(1.0-eta);
     62:   dN[0][1] =  0.25*(1.0-eta);
     63:   dN[0][2] =  0.25*(1.0+eta);
     64:   dN[0][3] = -0.25*(1.0+eta);

     66:   dN[1][0] = -0.25*(1.0-xi);
     67:   dN[1][1] = -0.25*(1.0+xi);
     68:   dN[1][2] =  0.25*(1.0+xi);
     69:   dN[1][3] =  0.25*(1.0-xi);
     70: }

     72: static void EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL],PetscScalar dNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
     73: {
     74:   PetscScalar J00,J01,J10,J11,J;
     75:   PetscScalar iJ00,iJ01,iJ10,iJ11;
     76:   PetscInt    i;

     78:   J00 = J01 = J10 = J11 = 0.0;
     79:   for (i = 0; i < NODES_PER_EL; i++) {
     80:     PetscScalar cx = coords[2*i];
     81:     PetscScalar cy = coords[2*i+1];

     83:     J00 += dN[0][i]*cx;      /* J_xx = dx/dxi */
     84:     J01 += dN[0][i]*cy;      /* J_xy = dy/dxi */
     85:     J10 += dN[1][i]*cx;      /* J_yx = dx/deta */
     86:     J11 += dN[1][i]*cy;      /* J_yy = dy/deta */
     87:   }
     88:   J = (J00*J11)-(J01*J10);

     90:   iJ00 =  J11/J;
     91:   iJ01 = -J01/J;
     92:   iJ10 = -J10/J;
     93:   iJ11 =  J00/J;

     95:   for (i = 0; i < NODES_PER_EL; i++) {
     96:     dNx[0][i] = dN[0][i]*iJ00+dN[1][i]*iJ01;
     97:     dNx[1][i] = dN[0][i]*iJ10+dN[1][i]*iJ11;
     98:   }

    100:   if (det_J) *det_J = J;
    101: }

    103: static void CreateGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
    104: {
    105:   *ngp         = 4;
    106:   gp_xi[0][0]  = -0.57735026919; gp_xi[0][1] = -0.57735026919;
    107:   gp_xi[1][0]  = -0.57735026919; gp_xi[1][1] =  0.57735026919;
    108:   gp_xi[2][0]  =  0.57735026919; gp_xi[2][1] =  0.57735026919;
    109:   gp_xi[3][0]  =  0.57735026919; gp_xi[3][1] = -0.57735026919;
    110:   gp_weight[0] = 1.0;
    111:   gp_weight[1] = 1.0;
    112:   gp_weight[2] = 1.0;
    113:   gp_weight[3] = 1.0;
    114: }

    116: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[],PetscInt s_u[],PetscInt s_p[])
    117: {
    118:   PetscInt i;
    120:   for (i=0; i<NODES_PER_EL; i++) {
    121:     /* velocity */
    122:     s_u[NSD*i+0] = 3*element[i];
    123:     s_u[NSD*i+1] = 3*element[i]+1;
    124:     /* pressure */
    125:     s_p[i] = 3*element[i]+2;
    126:   }
    127:   return(0);
    128: }

    130: static PetscInt map_wIwDI_uJuDJ(PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
    131: {
    132:   PetscInt ij,r,c,nc;

    134:   nc = u_NPE*u_dof;
    135:   r = w_dof*wi+wd;
    136:   c = u_dof*ui+ud;
    137:   ij = r*nc+c;
    138:   return(ij);
    139: }

    141: static void BForm_DivT(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
    142: {
    143:   PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
    144:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    145:   PetscScalar J_p,tildeD[3];
    146:   PetscScalar B[3][U_DOFS*NODES_PER_EL];
    147:   PetscInt    p,i,j,k,ngp;

    149:   /* define quadrature rule */
    150:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

    152:   /* evaluate bilinear form */
    153:   for (p = 0; p < ngp; p++) {
    154:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
    155:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);

    157:     for (i = 0; i < NODES_PER_EL; i++) {
    158:       PetscScalar d_dx_i = GNx_p[0][i];
    159:       PetscScalar d_dy_i = GNx_p[1][i];

    161:       B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
    162:       B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
    163:       B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
    164:     }

    166:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
    167:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
    168:     tildeD[2] =       gp_weight[p]*J_p*eta[p];

    170:     /* form Bt tildeD B */
    171:     /*
    172:     Ke_ij = Bt_ik . D_kl . B_lj
    173:           = B_ki . D_kl . B_lj
    174:           = B_ki . D_kk . B_kj
    175:     */
    176:     for (i = 0; i < 8; i++) {
    177:       for (j = 0; j < 8; j++) {
    178:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
    179:           Ke[i+8*j] += B[k][i]*tildeD[k]*B[k][j];
    180:         }
    181:       }
    182:     }
    183:   }
    184: }

    186: static void BForm_Grad(PetscScalar Ke[],PetscScalar coords[])
    187: {
    188:   PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
    189:   PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    190:   PetscScalar J_p,fac;
    191:   PetscInt    p,i,j,di,ngp;

    193:   /* define quadrature rule */
    194:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

    196:   /* evaluate bilinear form */
    197:   for (p = 0; p < ngp; p++) {
    198:     EvaluateBasis_Q1(gp_xi[p],Ni_p);
    199:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
    200:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
    201:     fac = gp_weight[p]*J_p;

    203:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
    204:       for (di = 0; di < NSD; di++) { /* u dofs */
    205:         for (j = 0; j < 4; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
    206:           PetscInt IJ;
    207:           IJ = map_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);

    209:           Ke[IJ] -= GNx_p[di][i]*Ni_p[j]*fac;
    210:         }
    211:       }
    212:     }
    213:   }
    214: }

    216: static void BForm_Div(PetscScalar De[],PetscScalar coords[])
    217: {
    218:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
    219:   PetscInt    i,j,nr_g,nc_g;

    221:   PetscMemzero(Ge,sizeof(Ge));
    222:   BForm_Grad(Ge,coords);

    224:   nr_g = U_DOFS*NODES_PER_EL;
    225:   nc_g = P_DOFS*NODES_PER_EL;

    227:   for (i = 0; i < nr_g; i++) {
    228:     for (j = 0; j < nc_g; j++) {
    229:       De[nr_g*j+i] = Ge[nc_g*i+j];
    230:     }
    231:   }
    232: }

    234: static void BForm_Stabilisation(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
    235: {
    236:   PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
    237:   PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    238:   PetscScalar J_p,fac,eta_avg;
    239:   PetscInt    p,i,j,ngp;

    241:   /* define quadrature rule */
    242:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

    244:   /* evaluate bilinear form */
    245:   for (p = 0; p < ngp; p++) {
    246:     EvaluateBasis_Q1(gp_xi[p],Ni_p);
    247:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
    248:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
    249:     fac = gp_weight[p]*J_p;

    251:     for (i = 0; i < NODES_PER_EL; i++) {
    252:       for (j = 0; j < NODES_PER_EL; j++) {
    253:         Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.0625);
    254:       }
    255:     }
    256:   }

    258:   /* scale */
    259:   eta_avg = 0.0;
    260:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    261:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
    262:   fac     = 1.0/eta_avg;
    263:   for (i = 0; i < NODES_PER_EL; i++) {
    264:     for (j = 0; j < NODES_PER_EL; j++) {
    265:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
    266:     }
    267:   }
    268: }

    270: static void BForm_ScaledMassMatrix(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
    271: {
    272:   PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
    273:   PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    274:   PetscScalar J_p,fac,eta_avg;
    275:   PetscInt    p,i,j,ngp;

    277:   /* define quadrature rule */
    278:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

    280:   /* evaluate bilinear form */
    281:   for (p = 0; p < ngp; p++) {
    282:     EvaluateBasis_Q1(gp_xi[p],Ni_p);
    283:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
    284:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
    285:     fac = gp_weight[p]*J_p;

    287:     for (i = 0; i < NODES_PER_EL; i++) {
    288:       for (j = 0; j < NODES_PER_EL; j++) {
    289:         Ke[NODES_PER_EL*i+j] -= fac*Ni_p[i]*Ni_p[j];
    290:       }
    291:     }
    292:   }

    294:   /* scale */
    295:   eta_avg = 0.0;
    296:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
    297:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
    298:   fac     = 1.0/eta_avg;
    299:   for (i = 0; i < NODES_PER_EL; i++) {
    300:     for (j = 0; j < NODES_PER_EL; j++) {
    301:       Ke[NODES_PER_EL*i+j] *= fac;
    302:     }
    303:   }
    304: }

    306: static void LForm_MomentumRHS(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
    307: {
    308:   PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
    309:   PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
    310:   PetscScalar J_p,fac;
    311:   PetscInt    p,i,ngp;

    313:   /* define quadrature rule */
    314:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

    316:   /* evaluate linear form */
    317:   for (p = 0; p < ngp; p++) {
    318:     EvaluateBasis_Q1(gp_xi[p],Ni_p);
    319:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
    320:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
    321:     fac = gp_weight[p]*J_p;

    323:     for (i = 0; i < NODES_PER_EL; i++) {
    324:       Fe[NSD*i]    = 0.0;
    325:       Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
    326:     }
    327:   }
    328: }

    330: static PetscErrorCode GetElementCoords(const PetscScalar _coords[],const PetscInt e2n[],PetscScalar el_coords[])
    331: {
    332:   PetscInt i,d;
    334:   /* get coords for the element */
    335:   for (i=0; i<4; i++) {
    336:     for (d=0; d<NSD; d++) {
    337:       el_coords[NSD*i+d] = _coords[NSD * e2n[i] + d];
    338:     }
    339:   }
    340:   return(0);
    341: }

    343: static PetscErrorCode AssembleStokes_A(Mat A,DM stokes_da,DM quadrature)
    344: {
    345:   DM                     cda;
    346:   Vec                    coords;
    347:   const PetscScalar      *_coords;
    348:   PetscInt               u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    349:   PetscInt               p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
    350:   PetscInt               nel,npe,eidx;
    351:   const PetscInt         *element_list;
    352:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
    353:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
    354:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
    355:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
    356:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    357:   PetscScalar            *q_eta,*prop_eta;
    358:   PetscErrorCode         ierr;

    361:   MatZeroEntries(A);
    362:   /* setup for coords */
    363:   DMGetCoordinateDM(stokes_da,&cda);
    364:   DMGetCoordinatesLocal(stokes_da,&coords);
    365:   VecGetArrayRead(coords,&_coords);

    367:   /* setup for coefficients */
    368:   DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);

    370:   DMDAGetElements(stokes_da,&nel,&npe,&element_list);
    371:   for (eidx = 0; eidx < nel; eidx++) {
    372:     const PetscInt *element = &element_list[npe*eidx];

    374:     /* get coords for the element */
    375:     GetElementCoords(_coords,element,el_coords);

    377:     /* get coefficients for the element */
    378:     prop_eta = &q_eta[GAUSS_POINTS * eidx];

    380:     /* initialise element stiffness matrix */
    381:     PetscMemzero(Ae,sizeof(Ae));
    382:     PetscMemzero(Ge,sizeof(Ge));
    383:     PetscMemzero(De,sizeof(De));
    384:     PetscMemzero(Ce,sizeof(Ce));

    386:     /* form element stiffness matrix */
    387:     BForm_DivT(Ae,el_coords,prop_eta);
    388:     BForm_Grad(Ge,el_coords);
    389:     BForm_Div(De,el_coords);
    390:     BForm_Stabilisation(Ce,el_coords,prop_eta);

    392:     /* insert element matrix into global matrix */
    393:     DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
    394:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
    395:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
    396:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
    397:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
    398:   }
    399:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    400:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    402:   DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
    403:   VecRestoreArrayRead(coords,&_coords);
    404:   return(0);
    405: }

    407: static PetscErrorCode AssembleStokes_PC(Mat A,DM stokes_da,DM quadrature)
    408: {
    409:   DM                     cda;
    410:   Vec                    coords;
    411:   const PetscScalar      *_coords;
    412:   PetscInt               u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    413:   PetscInt               p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
    414:   PetscInt               nel,npe,eidx;
    415:   const PetscInt         *element_list;
    416:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
    417:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
    418:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
    419:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
    420:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    421:   PetscScalar            *q_eta,*prop_eta;
    422:   PetscErrorCode         ierr;

    425:   MatZeroEntries(A);
    426:   /* setup for coords */
    427:   DMGetCoordinateDM(stokes_da,&cda);
    428:   DMGetCoordinatesLocal(stokes_da,&coords);
    429:   VecGetArrayRead(coords,&_coords);

    431:   /* setup for coefficients */
    432:   DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);

    434:   DMDAGetElements(stokes_da,&nel,&npe,&element_list);
    435:   for (eidx = 0; eidx < nel; eidx++) {
    436:     const PetscInt *element = &element_list[npe*eidx];

    438:     /* get coords for the element */
    439:     GetElementCoords(_coords,element,el_coords);

    441:     /* get coefficients for the element */
    442:     prop_eta = &q_eta[GAUSS_POINTS * eidx];

    444:     /* initialise element stiffness matrix */
    445:     PetscMemzero(Ae,sizeof(Ae));
    446:     PetscMemzero(Ge,sizeof(Ge));
    447:     PetscMemzero(De,sizeof(De));
    448:     PetscMemzero(Ce,sizeof(Ce));

    450:     /* form element stiffness matrix */
    451:     BForm_DivT(Ae,el_coords,prop_eta);
    452:     BForm_Grad(Ge,el_coords);
    453:     BForm_ScaledMassMatrix(Ce,el_coords,prop_eta);

    455:     /* insert element matrix into global matrix */
    456:     DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
    457:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
    458:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
    459:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
    460:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
    461:   }
    462:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    463:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    465:   DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
    466:   VecRestoreArrayRead(coords,&_coords);

    468:   return(0);
    469: }

    471: static PetscErrorCode AssembleStokes_RHS(Vec F,DM stokes_da,DM quadrature)
    472: {
    473:   DM                     cda;
    474:   Vec                    coords;
    475:   const PetscScalar      *_coords;
    476:   PetscInt               u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    477:   PetscInt               p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
    478:   PetscInt               nel,npe,eidx,i;
    479:   const PetscInt         *element_list;
    480:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
    481:   PetscScalar            He[NODES_PER_EL*P_DOFS];
    482:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    483:   PetscScalar            *q_rhs,*prop_fy;
    484:   Vec                    local_F;
    485:   PetscScalar            *LA_F;
    486:   PetscErrorCode         ierr;

    489:   VecZeroEntries(F);
    490:   /* setup for coords */
    491:   DMGetCoordinateDM(stokes_da,&cda);
    492:   DMGetCoordinatesLocal(stokes_da,&coords);
    493:   VecGetArrayRead(coords,&_coords);

    495:   /* setup for coefficients */
    496:   DMSwarmGetField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);

    498:   /* get acces to the vector */
    499:   DMGetLocalVector(stokes_da,&local_F);
    500:   VecZeroEntries(local_F);
    501:   VecGetArray(local_F,&LA_F);

    503:   DMDAGetElements(stokes_da,&nel,&npe,&element_list);
    504:   for (eidx = 0; eidx < nel; eidx++) {
    505:     const PetscInt *element = &element_list[npe*eidx];

    507:     /* get coords for the element */
    508:     GetElementCoords(_coords,element,el_coords);

    510:     /* get coefficients for the element */
    511:     prop_fy = &q_rhs[GAUSS_POINTS * eidx];

    513:     /* initialise element stiffness matrix */
    514:     PetscMemzero(Fe,sizeof(Fe));
    515:     PetscMemzero(He,sizeof(He));

    517:     /* form element stiffness matrix */
    518:     LForm_MomentumRHS(Fe,el_coords,NULL,prop_fy);

    520:     /* insert element matrix into global matrix */
    521:     DMDAGetElementEqnums_up(element,u_eqn,p_eqn);

    523:     for (i=0; i<NODES_PER_EL*U_DOFS; i++) {
    524:       LA_F[ u_eqn[i] ] += Fe[i];
    525:     }
    526:   }
    527:   DMSwarmRestoreField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
    528:   VecRestoreArrayRead(coords,&_coords);

    530:   VecRestoreArray(local_F,&LA_F);
    531:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
    532:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
    533:   DMRestoreLocalVector(stokes_da,&local_F);

    535:   return(0);
    536: }

    538: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm,DM dmc,PetscInt e,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization)
    539: {
    540:   PetscErrorCode    ierr;
    541:   PetscInt          dim,nel,npe,q,k,d,ncurr;
    542:   const PetscInt    *element_list;
    543:   Vec               coor;
    544:   const PetscScalar *_coor;
    545:   PetscReal         **basis,*elcoor,*xp;
    546:   PetscReal         *swarm_coor;
    547:   PetscInt          *swarm_cellid;

    550:   DMGetDimension(dm,&dim);
    551:   DMDAGetElements(dmc,&nel,&npe,&element_list);

    553:   PetscMalloc1(dim*npoints,&xp);
    554:   PetscMalloc1(dim*npe,&elcoor);
    555:   PetscMalloc1(npoints,&basis);
    556:   for (q=0; q<npoints; q++) {
    557:     PetscMalloc1(npe,&basis[q]);

    559:     switch (dim) {
    560:       case 1:
    561:         basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
    562:         basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
    563:         break;
    564:       case 2:
    565:         basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
    566:         basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
    567:         basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
    568:         basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
    569:         break;

    571:       case 3:
    572:         basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
    573:         basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
    574:         basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
    575:         basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
    576:         basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
    577:         basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
    578:         basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
    579:         basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
    580:         break;
    581:     }
    582:   }

    584:   DMGetCoordinatesLocal(dmc,&coor);
    585:   VecGetArrayRead(coor,&_coor);
    586:   /* compute and store the coordinates for the new points */
    587:   {
    588:     const PetscInt *element = &element_list[npe*e];

    590:     for (k=0; k<npe; k++) {
    591:       for (d=0; d<dim; d++) {
    592:         elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
    593:       }
    594:     }
    595:     for (q=0; q<npoints; q++) {
    596:       for (d=0; d<dim; d++) {
    597:         xp[dim*q+d] = 0.0;
    598:       }
    599:       for (k=0; k<npe; k++) {
    600:         for (d=0; d<dim; d++) {
    601:           xp[dim*q+d] += basis[q][k] * elcoor[dim*k+d];
    602:         }
    603:       }
    604:     }
    605:   }
    606:   VecRestoreArrayRead(coor,&_coor);
    607:   DMDARestoreElements(dmc,&nel,&npe,&element_list);

    609:   DMSwarmGetLocalSize(dm,&ncurr);
    610:   DMSwarmAddNPoints(dm,npoints);

    612:   if (proximity_initialization) {
    613:     PetscInt  *nnlist;
    614:     PetscReal *coor_q,*coor_qn;
    615:     PetscInt  npoints_e,*plist_e;

    617:     DMSwarmSortGetPointsPerCell(dm,e,&npoints_e,&plist_e);

    619:     PetscMalloc1(npoints,&nnlist);
    620:     /* find nearest neighour points in this cell */
    621:     DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    622:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    623:     for (q=0; q<npoints; q++) {
    624:       PetscInt  qn,nearest_neighour = -1;
    625:       PetscReal sep,min_sep = PETSC_MAX_REAL;

    627:       coor_q = &xp[dim*q];
    628:       for (qn=0; qn<npoints_e; qn++) {
    629:         coor_qn = &swarm_coor[dim*plist_e[qn]];
    630:         sep = 0.0;
    631:         for (d=0; d<dim; d++) {
    632:           sep += (coor_q[d]-coor_qn[d])*(coor_q[d]-coor_qn[d]);
    633:         }
    634:         if (sep < min_sep) {
    635:           nearest_neighour = plist_e[qn];
    636:           min_sep = sep;
    637:         }
    638:       }
    639:       if (nearest_neighour == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cell %D is empty - cannot initalize using nearest neighbours",e);
    640:       nnlist[q] = nearest_neighour;
    641:     }
    642:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    643:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

    645:     /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
    646:     for (q=0; q<npoints; q++) {
    647:       DMSwarmCopyPoint(dm,nnlist[q],ncurr+q);
    648:     }
    649:     DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    650:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    651:     for (q=0; q<npoints; q++) {
    652:       /* set the coordinates */
    653:       for (d=0; d<dim; d++) {
    654:         swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
    655:       }
    656:       /* set the cell index */
    657:       swarm_cellid[ncurr+q] = e;
    658:     }
    659:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    660:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

    662:     PetscFree(plist_e);
    663:     PetscFree(nnlist);
    664:   } else {
    665:     DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    666:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    667:     for (q=0; q<npoints; q++) {
    668:       /* set the coordinates */
    669:       for (d=0; d<dim; d++) {
    670:         swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
    671:       }
    672:       /* set the cell index */
    673:       swarm_cellid[ncurr+q] = e;
    674:     }
    675:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    676:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    677:   }

    679:   PetscFree(xp);
    680:   PetscFree(elcoor);
    681:   for (q=0; q<npoints; q++) {
    682:     PetscFree(basis[q]);
    683:   }
    684:   PetscFree(basis);
    685:   return(0);
    686: }

    688: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp,DM dm_mpoint)
    689: {
    690:   PetscInt        _npe,_nel,e,nel;
    691:   const PetscInt  *element;
    692:   DM              dmc;
    693:   PetscQuadrature quadrature;
    694:   const PetscReal *xi;
    695:   PetscInt        npoints_q,cnt,cnt_g;
    696:   PetscErrorCode  ierr;

    699:   DMDAGetElements(dm_vp,&_nel,&_npe,&element);
    700:   nel = _nel;
    701:   DMDARestoreElements(dm_vp,&_nel,&_npe,&element);

    703:   PetscDTGaussTensorQuadrature(2,1,4,-1.0,1.0,&quadrature);
    704:   PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&xi,NULL);
    705:   DMSwarmGetCellDM(dm_mpoint,&dmc);

    707:   DMSwarmSortGetAccess(dm_mpoint);

    709:   cnt = 0;
    710:   for (e=0; e<nel; e++) {
    711:     PetscInt npoints_per_cell;

    713:     DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint,e,&npoints_per_cell);

    715:     if (npoints_per_cell < 12) {
    716:       DMSwarmPICInsertPointsCellwise(dm_mpoint,dm_vp,e,npoints_q,(PetscReal*)xi,PETSC_TRUE);
    717:       cnt++;
    718:     }
    719:   }
    720:   MPI_Allreduce(&cnt,&cnt_g,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
    721:   if (cnt_g > 0) {
    722:     PetscPrintf(PETSC_COMM_WORLD,".... ....pop cont: adjusted %D cells\n",cnt_g);
    723:   }

    725:   DMSwarmSortRestoreAccess(dm_mpoint);
    726:   PetscQuadratureDestroy(&quadrature);
    727:   return(0);
    728: }

    730: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp,Vec vp,PetscReal dt,DM dm_mpoint)
    731: {
    732:   PetscErrorCode    ierr;
    733:   Vec               vp_l,coor_l;
    734:   const PetscScalar *LA_vp;
    735:   PetscInt          i,p,e,npoints,nel,npe;
    736:   PetscInt          *mpfield_cell;
    737:   PetscReal         *mpfield_coor;
    738:   const PetscInt    *element_list;
    739:   const PetscInt    *element;
    740:   PetscScalar       xi_p[NSD],Ni[NODES_PER_EL];
    741:   const PetscScalar *LA_coor;
    742:   PetscScalar       dx[NSD];

    745:   DMGetCoordinatesLocal(dm_vp,&coor_l);
    746:   VecGetArrayRead(coor_l,&LA_coor);

    748:   DMGetLocalVector(dm_vp,&vp_l);
    749:   DMGlobalToLocalBegin(dm_vp,vp,INSERT_VALUES,vp_l);
    750:   DMGlobalToLocalEnd(dm_vp,vp,INSERT_VALUES,vp_l);
    751:   VecGetArrayRead(vp_l,&LA_vp);

    753:   DMDAGetElements(dm_vp,&nel,&npe,&element_list);
    754:   DMSwarmGetLocalSize(dm_mpoint,&npoints);
    755:   DMSwarmGetField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
    756:   DMSwarmGetField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
    757:   for (p=0; p<npoints; p++) {
    758:     PetscReal         *coor_p;
    759:     PetscScalar       vel_n[NSD*NODES_PER_EL],vel_p[NSD];
    760:     const PetscScalar *x0;
    761:     const PetscScalar *x2;

    763:     e       = mpfield_cell[p];
    764:     coor_p  = &mpfield_coor[NSD*p];
    765:     element = &element_list[NODES_PER_EL*e];

    767:     /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
    768:     x0 = &LA_coor[NSD*element[0]];
    769:     x2 = &LA_coor[NSD*element[2]];

    771:     dx[0] = x2[0] - x0[0];
    772:     dx[1] = x2[1] - x0[1];

    774:     xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
    775:     xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;
    776:     if (PetscRealPart(xi_p[0]) < -1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (xi) too small %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[0]),e);
    777:     if (PetscRealPart(xi_p[0]) >  1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (xi) too large %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[0]),e);
    778:     if (PetscRealPart(xi_p[1]) < -1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (eta) too small %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[1]),e);
    779:     if (PetscRealPart(xi_p[1]) >  1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (eta) too large %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[1]),e);

    781:     /* evaluate basis functions */
    782:     EvaluateBasis_Q1(xi_p,Ni);

    784:     /* get cell nodal velocities */
    785:     for (i=0; i<NODES_PER_EL; i++) {
    786:       PetscInt nid;

    788:       nid = element[i];
    789:       vel_n[NSD*i+0] = LA_vp[(NSD+1)*nid+0];
    790:       vel_n[NSD*i+1] = LA_vp[(NSD+1)*nid+1];
    791:     }

    793:     /* interpolate velocity */
    794:     vel_p[0] = vel_p[1] = 0.0;
    795:     for (i=0; i<NODES_PER_EL; i++) {
    796:       vel_p[0] += Ni[i] * vel_n[NSD*i+0];
    797:       vel_p[1] += Ni[i] * vel_n[NSD*i+1];
    798:     }

    800:     coor_p[0] += dt * PetscRealPart(vel_p[0]);
    801:     coor_p[1] += dt * PetscRealPart(vel_p[1]);
    802:   }

    804:   DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
    805:   DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
    806:   DMDARestoreElements(dm_vp,&nel,&npe,&element_list);
    807:   VecRestoreArrayRead(vp_l,&LA_vp);
    808:   DMRestoreLocalVector(dm_vp,&vp_l);
    809:   VecRestoreArrayRead(coor_l,&LA_coor);
    810:   return(0);
    811: }

    813: PetscErrorCode MaterialPoint_Interpolate(DM dm,Vec eta_v,Vec rho_v,DM dm_quadrature)
    814: {
    815:   Vec            eta_l,rho_l;
    816:   PetscScalar    *_eta_l,*_rho_l;
    817:   PetscInt       nqp,npe,nel;
    818:   PetscScalar    qp_xi[GAUSS_POINTS][NSD];
    819:   PetscScalar    qp_weight[GAUSS_POINTS];
    820:   PetscInt       q,k,e;
    821:   PetscScalar    Ni[GAUSS_POINTS][NODES_PER_EL];
    822:   const PetscInt *element_list;
    823:   PetscReal      *q_eta,*q_rhs;

    827:   /* define quadrature rule */
    828:   CreateGaussQuadrature(&nqp,qp_xi,qp_weight);
    829:   for (q=0; q<nqp; q++) {
    830:     EvaluateBasis_Q1(qp_xi[q],Ni[q]);
    831:   }

    833:   DMGetLocalVector(dm,&eta_l);
    834:   DMGetLocalVector(dm,&rho_l);

    836:   DMGlobalToLocalBegin(dm,eta_v,INSERT_VALUES,eta_l);
    837:   DMGlobalToLocalEnd(dm,eta_v,INSERT_VALUES,eta_l);
    838:   DMGlobalToLocalBegin(dm,rho_v,INSERT_VALUES,rho_l);
    839:   DMGlobalToLocalEnd(dm,rho_v,INSERT_VALUES,rho_l);

    841:   VecGetArray(eta_l,&_eta_l);
    842:   VecGetArray(rho_l,&_rho_l);

    844:   DMSwarmGetField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
    845:   DMSwarmGetField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);

    847:   DMDAGetElements(dm,&nel,&npe,&element_list);
    848:   for (e=0; e<nel; e++) {
    849:     PetscScalar    eta_field_e[NODES_PER_EL];
    850:     PetscScalar    rho_field_e[NODES_PER_EL];
    851:     const PetscInt *element = &element_list[4*e];

    853:     for (k=0; k<NODES_PER_EL; k++) {
    854:       eta_field_e[k] = _eta_l[ element[k] ];
    855:       rho_field_e[k] = _rho_l[ element[k] ];
    856:     }

    858:     for (q=0; q<nqp; q++) {
    859:       PetscScalar eta_q,rho_q;

    861:       eta_q = rho_q = 0.0;
    862:       for (k=0; k<NODES_PER_EL; k++) {
    863:         eta_q += Ni[q][k] * eta_field_e[k];
    864:         rho_q += Ni[q][k] * rho_field_e[k];
    865:       }

    867:       q_eta[nqp*e+q] = PetscRealPart(eta_q);
    868:       q_rhs[nqp*e+q] = PetscRealPart(rho_q);
    869:     }
    870:   }
    871:   DMDARestoreElements(dm,&nel,&npe,&element_list);

    873:   DMSwarmRestoreField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
    874:   DMSwarmRestoreField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);

    876:   VecRestoreArray(rho_l,&_rho_l);
    877:   VecRestoreArray(eta_l,&_eta_l);
    878:   DMRestoreLocalVector(dm,&rho_l);
    879:   DMRestoreLocalVector(dm,&eta_l);
    880:   return(0);
    881: }

    883: static PetscErrorCode SolveTimeDepStokes(PetscInt mx,PetscInt my)
    884: {
    885:   DM                     dm_stokes,dm_coeff;
    886:   PetscInt               u_dof,p_dof,dof,stencil_width;
    887:   Mat                    A,B;
    888:   PetscInt               nel_local;
    889:   Vec                    eta_v,rho_v;
    890:   Vec                    f,X;
    891:   KSP                    ksp;
    892:   PC                     pc;
    893:   char                   filename[PETSC_MAX_PATH_LEN];
    894:   DM                     dms_quadrature,dms_mpoint;
    895:   PetscInt               nel,npe,npoints;
    896:   const PetscInt         *element_list;
    897:   PetscInt               tk,nt,dump_freq;
    898:   PetscReal              dt,dt_max = 0.0;
    899:   PetscReal              vx[2],vy[2],max_v = 0.0,max_v_step,dh;
    900:   PetscErrorCode         ierr;
    901:   const char             *fieldnames[] = { "eta" , "rho" };
    902:   Vec                    *pfields;
    903:   PetscInt               ppcell = 1;
    904:   PetscReal              time,delta_eta = 1.0;
    905:   PetscBool              randomize_coords = PETSC_FALSE;
    906:   PetscReal              randomize_fac = 0.25;
    907:   PetscBool              no_view = PETSC_FALSE;
    908:   PetscBool              isbddc;

    911:   /*
    912:     Generate the DMDA for the velocity and pressure spaces.
    913:     We use Q1 elements for both fields.
    914:     The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
    915:     The number of nodes in each direction is mx+1, my+1
    916:   */
    917:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
    918:   p_dof         = P_DOFS; /* p - pressure */
    919:   dof           = u_dof + p_dof;
    920:   stencil_width = 1;
    921:   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,&dm_stokes);
    922:   DMDASetElementType(dm_stokes,DMDA_ELEMENT_Q1);
    923:   DMSetMatType(dm_stokes,MATAIJ);
    924:   DMSetFromOptions(dm_stokes);
    925:   DMSetUp(dm_stokes);
    926:   DMDASetFieldName(dm_stokes,0,"ux");
    927:   DMDASetFieldName(dm_stokes,1,"uy");
    928:   DMDASetFieldName(dm_stokes,2,"p");

    930:   /* unit box [0,0.9142] x [0,1] */
    931:   DMDASetUniformCoordinates(dm_stokes,0.0,0.9142,0.0,1.0,0.,0.);
    932:   dh = 1.0/((PetscReal)(mx));

    934:   /* Get local number of elements */
    935:   {
    936:     DMDAGetElements(dm_stokes,&nel,&npe,&element_list);

    938:     nel_local = nel;

    940:     DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
    941:   }

    943:   /* Create DMDA for representing scalar fields */
    944:   DMDACreateCompatibleDMDA(dm_stokes,1,&dm_coeff);

    946:   /* Create the swarm for storing quadrature point values */
    947:   DMCreate(PETSC_COMM_WORLD,&dms_quadrature);
    948:   DMSetType(dms_quadrature,DMSWARM);
    949:   DMSetDimension(dms_quadrature,2);

    951:   /* Register fields for viscosity and density on the quadrature points */
    952:   DMSwarmRegisterPetscDatatypeField(dms_quadrature,"eta_q",1,PETSC_REAL);
    953:   DMSwarmRegisterPetscDatatypeField(dms_quadrature,"rho_q",1,PETSC_REAL);
    954:   DMSwarmFinalizeFieldRegister(dms_quadrature);
    955:   DMSwarmSetLocalSizes(dms_quadrature,nel_local * GAUSS_POINTS,0);

    957:   /* Create the material point swarm */
    958:   DMCreate(PETSC_COMM_WORLD,&dms_mpoint);
    959:   DMSetType(dms_mpoint,DMSWARM);
    960:   DMSetDimension(dms_mpoint,2);

    962:   /* Configure the material point swarm to be of type Particle-In-Cell */
    963:   DMSwarmSetType(dms_mpoint,DMSWARM_PIC);

    965:   /*
    966:      Specify the DM to use for point location and projections
    967:      within the context of a PIC scheme
    968:   */
    969:   DMSwarmSetCellDM(dms_mpoint,dm_coeff);

    971:   /* Register fields for viscosity and density */
    972:   DMSwarmRegisterPetscDatatypeField(dms_mpoint,"eta",1,PETSC_REAL);
    973:   DMSwarmRegisterPetscDatatypeField(dms_mpoint,"rho",1,PETSC_REAL);
    974:   DMSwarmFinalizeFieldRegister(dms_mpoint);

    976:   PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);
    977:   DMSwarmSetLocalSizes(dms_mpoint,nel_local * ppcell,100);

    979:   /*
    980:     Layout the material points in space using the cell DM.
    981:     Particle coordinates are defined by cell wise using different methods.
    982:     - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
    983:                               corresponding to a Gauss quadrature rule with
    984:                               ppcell points in each direction.
    985:     - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
    986:                                 ppcell x ppcell quadralaterals defined within the
    987:                                 reference element.
    988:     - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
    989:                                     of each quadralateral obtained by sub-dividing
    990:                                     the reference element cell ppcell times.
    991:   */
    992:   DMSwarmInsertPointsUsingCellDM(dms_mpoint,DMSWARMPIC_LAYOUT_SUBDIVISION,ppcell);

    994:   /*
    995:     Defne a high resolution layer of material points across the material interface
    996:   */
    997:   {
    998:     PetscInt  npoints_dir_x[2];
    999:     PetscReal min[2],max[2];

    1001:     npoints_dir_x[0] = (PetscInt)(0.9142/(0.05*dh));
    1002:     npoints_dir_x[1] = (PetscInt)((0.25-0.15)/(0.05*dh));
    1003:     min[0] = 0.0;  max[0] = 0.9142;
    1004:     min[1] = 0.05; max[1] = 0.35;
    1005:     DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
    1006:   }

    1008:   /*
    1009:     Define a high resolution layer of material points near the surface of the domain
    1010:     to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
    1011:     when applied to buouyancy driven flow. The error in div(u) is O(h).
    1012:   */
    1013:   {
    1014:     PetscInt  npoints_dir_x[2];
    1015:     PetscReal min[2],max[2];

    1017:     npoints_dir_x[0] = (PetscInt)(0.9142/(0.25*dh));
    1018:     npoints_dir_x[1] = (PetscInt)(3.0*dh/(0.25*dh));
    1019:     min[0] = 0.0;          max[0] = 0.9142;
    1020:     min[1] = 1.0 - 3.0*dh; max[1] = 1.0-0.0001;
    1021:     DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
    1022:   }

    1024:   DMView(dms_mpoint,PETSC_VIEWER_STDOUT_WORLD);

    1026:   /* Define initial material properties on each particle in the material point swarm */
    1027:   PetscOptionsGetReal(NULL,NULL,"-delta_eta",&delta_eta,NULL);
    1028:   PetscOptionsGetBool(NULL,NULL,"-randomize_coords",&randomize_coords,NULL);
    1029:   PetscOptionsGetReal(NULL,NULL,"-randomize_fac",&randomize_fac,NULL);
    1030:   if (randomize_fac > 1.0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"The value of -randomize_fac should be <= 1.0");
    1031:   {
    1032:     PetscReal   *array_x,*array_e,*array_r;
    1033:     PetscInt    p;
    1034:     PetscRandom r;
    1035:     PetscMPIInt rank;

    1037:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

    1039:     PetscRandomCreate(PETSC_COMM_SELF,&r);
    1040:     PetscRandomSetInterval(r,-randomize_fac*dh,randomize_fac*dh);
    1041:     PetscRandomSetSeed(r,(unsigned long)rank);
    1042:     PetscRandomSeed(r);

    1044:     DMDAGetElements(dm_stokes,&nel,&npe,&element_list);

    1046:     /*
    1047:        Fetch the registered data from the material point DMSwarm.
    1048:        The fields "eta" and "rho" were registered by this example.
    1049:        The field identified by the the variable DMSwarmPICField_coor
    1050:        was registered by the DMSwarm implementation when the function
    1051:          DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
    1052:        was called. The returned array defines the coordinates of each
    1053:        material point in the point swarm.
    1054:     */
    1055:     DMSwarmGetField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);
    1056:     DMSwarmGetField(dms_mpoint,"eta",               NULL,NULL,(void**)&array_e);
    1057:     DMSwarmGetField(dms_mpoint,"rho",               NULL,NULL,(void**)&array_r);

    1059:     DMSwarmGetLocalSize(dms_mpoint,&npoints);
    1060:     for (p = 0; p < npoints; p++) {
    1061:       PetscReal x_p[2],rr[2];

    1063:       if (randomize_coords) {
    1064:         PetscRandomGetValueReal(r,&rr[0]);
    1065:         PetscRandomGetValueReal(r,&rr[1]);
    1066:         array_x[2*p + 0] += rr[0];
    1067:         array_x[2*p + 1] += rr[1];
    1068:       }

    1070:       /* Get the coordinates of point, p */
    1071:       x_p[0] = array_x[2*p + 0];
    1072:       x_p[1] = array_x[2*p + 1];

    1074:        if (x_p[1] < (0.2 + 0.02*PetscCosReal(PETSC_PI*x_p[0]/0.9142))) {
    1075:          /* Material properties below the interface */
    1076:          array_e[p] = 1.0 * (1.0/delta_eta);
    1077:          array_r[p] = 0.0;
    1078:        } else {
    1079:          /* Material properties above the interface */
    1080:          array_e[p] = 1.0;
    1081:          array_r[p] = 1.0;
    1082:        }
    1083:     }

    1085:     /*
    1086:        Restore the fetched data fields from the material point DMSwarm.
    1087:        Calling the Restore function invalidates the points array_r, array_e, array_x
    1088:        by setting them to NULL.
    1089:     */
    1090:     DMSwarmRestoreField(dms_mpoint,"rho",NULL,NULL,(void**)&array_r);
    1091:     DMSwarmRestoreField(dms_mpoint,"eta",NULL,NULL,(void**)&array_e);
    1092:     DMSwarmRestoreField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);

    1094:     DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
    1095:     PetscRandomDestroy(&r);
    1096:   }

    1098:   /*
    1099:      If the particle coordinates where randomly shifted, they may have crossed into another
    1100:      element, or into another sub-domain. To account for this we call the Migrate function.
    1101:   */
    1102:   if (randomize_coords) {
    1103:     DMSwarmMigrate(dms_mpoint,PETSC_TRUE);
    1104:   }

    1106:   PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
    1107:   if (!no_view) {
    1108:     DMSwarmViewXDMF(dms_mpoint,"ic_coeff_dms.xmf");
    1109:   }

    1111:   /* project the swarm properties */
    1112:   DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_FALSE);
    1113:   eta_v = pfields[0];
    1114:   rho_v = pfields[1];
    1115:   PetscObjectSetName((PetscObject)eta_v,"eta");
    1116:   PetscObjectSetName((PetscObject)rho_v,"rho");
    1117:   MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);

    1119:   /* view projected coefficients eta and rho */
    1120:   if (!no_view) {
    1121:     PetscViewer viewer;

    1123:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
    1124:     PetscViewerSetType(viewer,PETSCVIEWERVTK);
    1125:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
    1126:     PetscViewerFileSetName(viewer,"ic_coeff_dmda.vts");
    1127:     VecView(eta_v,viewer);
    1128:     VecView(rho_v,viewer);
    1129:     PetscViewerDestroy(&viewer);
    1130:   }

    1132:   DMCreateMatrix(dm_stokes,&A);
    1133:   DMCreateMatrix(dm_stokes,&B);
    1134:   DMCreateGlobalVector(dm_stokes,&f);
    1135:   DMCreateGlobalVector(dm_stokes,&X);

    1137:   AssembleStokes_A(A,dm_stokes,dms_quadrature);
    1138:   AssembleStokes_PC(B,dm_stokes,dms_quadrature);
    1139:   AssembleStokes_RHS(f,dm_stokes,dms_quadrature);

    1141:   DMDAApplyBoundaryConditions(dm_stokes,A,f);
    1142:   DMDAApplyBoundaryConditions(dm_stokes,B,NULL);

    1144:   KSPCreate(PETSC_COMM_WORLD,&ksp);
    1145:   KSPSetOptionsPrefix(ksp,"stokes_");
    1146:   KSPSetDM(ksp,dm_stokes);
    1147:   KSPSetDMActive(ksp,PETSC_FALSE);
    1148:   KSPSetOperators(ksp,A,B);
    1149:   KSPSetFromOptions(ksp);
    1150:   KSPGetPC(ksp,&pc);
    1151:   PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);
    1152:   if (isbddc) {
    1153:     KSPSetOperators(ksp,A,A);
    1154:   }

    1156:   /* Define u-v-p indices for fieldsplit */
    1157:   {
    1158:     PC             pc;
    1159:     const PetscInt ufields[] = {0,1},pfields[1] = {2};

    1161:     KSPGetPC(ksp,&pc);
    1162:     PCFieldSplitSetBlockSize(pc,3);
    1163:     PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
    1164:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
    1165:   }

    1167:   /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
    1168:   {
    1169:     PC        pc,pc_u;
    1170:     KSP       *sub_ksp,ksp_u;
    1171:     PetscInt  nsplits;
    1172:     DM        dm_u;
    1173:     PetscBool is_pcfs;

    1175:     KSPGetPC(ksp,&pc);

    1177:     is_pcfs = PETSC_FALSE;
    1178:     PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);

    1180:     if (is_pcfs) {
    1181:       KSPSetUp(ksp);
    1182:       KSPGetPC(ksp,&pc);
    1183:       PCFieldSplitGetSubKSP(pc,&nsplits,&sub_ksp);
    1184:       ksp_u = sub_ksp[0];
    1185:       PetscFree(sub_ksp);

    1187:       if (nsplits == 2) {
    1188:         DMDACreateCompatibleDMDA(dm_stokes,2,&dm_u);

    1190:         KSPSetDM(ksp_u,dm_u);
    1191:         KSPSetDMActive(ksp_u,PETSC_FALSE);
    1192:         DMDestroy(&dm_u);

    1194:         /* enforce galerkin coarse grids be used */
    1195:         KSPGetPC(ksp_u,&pc_u);
    1196:         PCMGSetGalerkin(pc_u,PC_MG_GALERKIN_PMAT);
    1197:       }
    1198:     }
    1199:   }

    1201:   dump_freq = 10;
    1202:   PetscOptionsGetInt(NULL,NULL,"-dump_freq",&dump_freq,NULL);
    1203:   nt = 10;
    1204:   PetscOptionsGetInt(NULL,NULL,"-nt",&nt,NULL);
    1205:   time = 0.0;
    1206:   for (tk=1; tk<=nt; tk++) {

    1208:     PetscPrintf(PETSC_COMM_WORLD,".... assemble\n");
    1209:     AssembleStokes_A(A,dm_stokes,dms_quadrature);
    1210:     AssembleStokes_PC(B,dm_stokes,dms_quadrature);
    1211:     AssembleStokes_RHS(f,dm_stokes,dms_quadrature);

    1213:     PetscPrintf(PETSC_COMM_WORLD,".... bc imposition\n");
    1214:     DMDAApplyBoundaryConditions(dm_stokes,A,f);
    1215:     DMDAApplyBoundaryConditions(dm_stokes,B,NULL);

    1217:     PetscPrintf(PETSC_COMM_WORLD,".... solve\n");
    1218:     KSPSetOperators(ksp,A, isbddc ? A : B);
    1219:     KSPSolve(ksp,f,X);

    1221:     VecStrideMax(X,0,NULL,&vx[1]);
    1222:     VecStrideMax(X,1,NULL,&vy[1]);
    1223:     VecStrideMin(X,0,NULL,&vx[0]);
    1224:     VecStrideMin(X,1,NULL,&vy[0]);

    1226:     max_v_step = PetscMax(vx[0],vx[1]);
    1227:     max_v_step = PetscMax(max_v_step,vy[0]);
    1228:     max_v_step = PetscMax(max_v_step,vy[1]);
    1229:     max_v = PetscMax(max_v,max_v_step);

    1231:     dt_max = 2.0;
    1232:     dt = 0.5 * (dh / max_v_step);
    1233:     PetscPrintf(PETSC_COMM_WORLD,".... max v %1.4e , dt %1.4e : [total] max v %1.4e , dt_max %1.4e\n",(double)max_v_step,(double)dt,(double)max_v,(double)dt_max);
    1234:     dt = PetscMin(dt_max,dt);

    1236:     /* advect */
    1237:     PetscPrintf(PETSC_COMM_WORLD,".... advect\n");
    1238:     MaterialPoint_AdvectRK1(dm_stokes,X,dt,dms_mpoint);

    1240:     /* migrate */
    1241:     PetscPrintf(PETSC_COMM_WORLD,".... migrate\n");
    1242:     DMSwarmMigrate(dms_mpoint,PETSC_TRUE);

    1244:     /* update cell population */
    1245:     PetscPrintf(PETSC_COMM_WORLD,".... populate cells\n");
    1246:     MaterialPoint_PopulateCell(dm_stokes,dms_mpoint);

    1248:     /* update coefficients on quadrature points */
    1249:     PetscPrintf(PETSC_COMM_WORLD,".... project\n");
    1250:     DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_TRUE);
    1251:     eta_v = pfields[0];
    1252:     rho_v = pfields[1];
    1253:     PetscPrintf(PETSC_COMM_WORLD,".... interp\n");
    1254:     MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);

    1256:     if (tk%dump_freq == 0) {
    1257:       PetscViewer viewer;

    1259:       PetscPrintf(PETSC_COMM_WORLD,".... write XDMF, VTS\n");
    1260:       PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_coeff_dms.xmf",tk);
    1261:       DMSwarmViewXDMF(dms_mpoint,filename);

    1263:       PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_vp_dm.vts",tk);
    1264:       PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
    1265:       PetscViewerSetType(viewer,PETSCVIEWERVTK);
    1266:       PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
    1267:       PetscViewerFileSetName(viewer,filename);
    1268:       VecView(X,viewer);
    1269:       PetscViewerDestroy(&viewer);
    1270:     }
    1271:     time += dt;
    1272:     PetscPrintf(PETSC_COMM_WORLD,"step %D : time %1.2e \n",tk,(double)time);
    1273:   }

    1275:   KSPDestroy(&ksp);
    1276:   VecDestroy(&X);
    1277:   VecDestroy(&f);
    1278:   MatDestroy(&A);
    1279:   MatDestroy(&B);
    1280:   VecDestroy(&eta_v);
    1281:   VecDestroy(&rho_v);
    1282:   PetscFree(pfields);

    1284:   DMDestroy(&dms_mpoint);
    1285:   DMDestroy(&dms_quadrature);
    1286:   DMDestroy(&dm_coeff);
    1287:   DMDestroy(&dm_stokes);
    1288:   return(0);
    1289: }

    1291: /*
    1292:  <sequential run>
    1293:  ./ex70 -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 lu  -mx 80 -my 80  -stokes_ksp_converged_reason  -dump_freq 25  -stokes_ksp_rtol 1.0e-8 -build_twosided allreduce  -ppcell 2 -nt 4000 -delta_eta 1.0 -randomize_coords
    1294: */
    1295: int main(int argc,char **args)
    1296: {
    1298:   PetscInt       mx,my;
    1299:   PetscBool      set = PETSC_FALSE;

    1301:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
    1302:   mx = my = 10;
    1303:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
    1304:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
    1305:   PetscOptionsGetInt(NULL,NULL,"-mxy",&mx,&set);
    1306:   if (set) {
    1307:     my = mx;
    1308:   }
    1309:   SolveTimeDepStokes(mx,my);
    1310:   PetscFinalize();
    1311:   return ierr;
    1312: }

    1314: /* -------------------------- helpers for boundary conditions -------------------------------- */
    1315: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
    1316: {
    1317:   DM                     cda;
    1318:   Vec                    coords;
    1319:   PetscInt               si,sj,nx,ny,i,j;
    1320:   PetscInt               M,N;
    1321:   DMDACoor2d             **_coords;
    1322:   const PetscInt         *g_idx;
    1323:   PetscInt               *bc_global_ids;
    1324:   PetscScalar            *bc_vals;
    1325:   PetscInt               nbcs;
    1326:   PetscInt               n_dofs;
    1327:   PetscErrorCode         ierr;
    1328:   ISLocalToGlobalMapping ltogm;

    1331:   DMGetLocalToGlobalMapping(da,&ltogm);
    1332:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1334:   DMGetCoordinateDM(da,&cda);
    1335:   DMGetCoordinatesLocal(da,&coords);
    1336:   DMDAVecGetArray(cda,coords,&_coords);
    1337:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1338:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1340:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
    1341:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

    1346:   i = nx-1;
    1347:   for (j = 0; j < ny; j++) {
    1348:     PetscInt local_id;

    1350:     local_id = i+j*nx;

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

    1354:     bc_vals[j] =  0.0;
    1355:   }
    1356:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1357:   nbcs = 0;
    1358:   if ((si+nx) == (M)) nbcs = ny;

    1360:   if (b) {
    1361:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1362:     VecAssemblyBegin(b);
    1363:     VecAssemblyEnd(b);
    1364:   }
    1365:   if (A) {
    1366:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1367:   }

    1369:   PetscFree(bc_vals);
    1370:   PetscFree(bc_global_ids);

    1372:   DMDAVecRestoreArray(cda,coords,&_coords);
    1373:   return(0);
    1374: }

    1376: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
    1377: {
    1378:   DM                     cda;
    1379:   Vec                    coords;
    1380:   PetscInt               si,sj,nx,ny,i,j;
    1381:   PetscInt               M,N;
    1382:   DMDACoor2d             **_coords;
    1383:   const PetscInt         *g_idx;
    1384:   PetscInt               *bc_global_ids;
    1385:   PetscScalar            *bc_vals;
    1386:   PetscInt               nbcs;
    1387:   PetscInt               n_dofs;
    1388:   PetscErrorCode         ierr;
    1389:   ISLocalToGlobalMapping ltogm;

    1392:   DMGetLocalToGlobalMapping(da,&ltogm);
    1393:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1395:   DMGetCoordinateDM(da,&cda);
    1396:   DMGetCoordinatesLocal(da,&coords);
    1397:   DMDAVecGetArray(cda,coords,&_coords);
    1398:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1399:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1401:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
    1402:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

    1407:   i = 0;
    1408:   for (j = 0; j < ny; j++) {
    1409:     PetscInt local_id;

    1411:     local_id = i+j*nx;

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

    1415:     bc_vals[j] =  0.0;
    1416:   }
    1417:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1418:   nbcs = 0;
    1419:   if (si == 0) nbcs = ny;

    1421:   if (b) {
    1422:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1423:     VecAssemblyBegin(b);
    1424:     VecAssemblyEnd(b);
    1425:   }

    1427:   if (A) {
    1428:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1429:   }

    1431:   PetscFree(bc_vals);
    1432:   PetscFree(bc_global_ids);

    1434:   DMDAVecRestoreArray(cda,coords,&_coords);
    1435:   return(0);
    1436: }

    1438: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
    1439: {
    1440:   DM                     cda;
    1441:   Vec                    coords;
    1442:   PetscInt               si,sj,nx,ny,i,j;
    1443:   PetscInt               M,N;
    1444:   DMDACoor2d             **_coords;
    1445:   const PetscInt         *g_idx;
    1446:   PetscInt               *bc_global_ids;
    1447:   PetscScalar            *bc_vals;
    1448:   PetscInt               nbcs;
    1449:   PetscInt               n_dofs;
    1450:   PetscErrorCode         ierr;
    1451:   ISLocalToGlobalMapping ltogm;

    1454:   DMGetLocalToGlobalMapping(da,&ltogm);
    1455:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1457:   DMGetCoordinateDM(da,&cda);
    1458:   DMGetCoordinatesLocal(da,&coords);
    1459:   DMDAVecGetArray(cda,coords,&_coords);
    1460:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1461:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1463:   PetscMalloc1(nx,&bc_global_ids);
    1464:   PetscMalloc1(nx,&bc_vals);

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

    1469:   j = ny-1;
    1470:   for (i = 0; i < nx; i++) {
    1471:     PetscInt local_id;

    1473:     local_id = i+j*nx;

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

    1477:     bc_vals[i] =  0.0;
    1478:   }
    1479:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1480:   nbcs = 0;
    1481:   if ((sj+ny) == (N)) nbcs = nx;

    1483:   if (b) {
    1484:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1485:     VecAssemblyBegin(b);
    1486:     VecAssemblyEnd(b);
    1487:   }
    1488:   if (A) {
    1489:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1490:   }

    1492:   PetscFree(bc_vals);
    1493:   PetscFree(bc_global_ids);

    1495:   DMDAVecRestoreArray(cda,coords,&_coords);
    1496:   return(0);
    1497: }

    1499: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
    1500: {
    1501:   DM                     cda;
    1502:   Vec                    coords;
    1503:   PetscInt               si,sj,nx,ny,i,j;
    1504:   PetscInt               M,N;
    1505:   DMDACoor2d             **_coords;
    1506:   const PetscInt         *g_idx;
    1507:   PetscInt               *bc_global_ids;
    1508:   PetscScalar            *bc_vals;
    1509:   PetscInt               nbcs;
    1510:   PetscInt               n_dofs;
    1511:   PetscErrorCode         ierr;
    1512:   ISLocalToGlobalMapping ltogm;

    1515:   DMGetLocalToGlobalMapping(da,&ltogm);
    1516:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1518:   DMGetCoordinateDM(da,&cda);
    1519:   DMGetCoordinatesLocal(da,&coords);
    1520:   DMDAVecGetArray(cda,coords,&_coords);
    1521:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1522:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1524:   PetscMalloc1(nx,&bc_global_ids);
    1525:   PetscMalloc1(nx,&bc_vals);

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

    1530:   j = 0;
    1531:   for (i = 0; i < nx; i++) {
    1532:     PetscInt local_id;

    1534:     local_id = i+j*nx;

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

    1538:     bc_vals[i] =  0.0;
    1539:   }
    1540:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1541:   nbcs = 0;
    1542:   if (sj == 0) nbcs = nx;

    1544:   if (b) {
    1545:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1546:     VecAssemblyBegin(b);
    1547:     VecAssemblyEnd(b);
    1548:   }
    1549:   if (A) {
    1550:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
    1551:   }

    1553:   PetscFree(bc_vals);
    1554:   PetscFree(bc_global_ids);

    1556:   DMDAVecRestoreArray(cda,coords,&_coords);
    1557:   return(0);
    1558: }

    1560: /*
    1561:  Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
    1562:  Impose no slip boundray conditions on the top/bottom faces:   u_i n_i = 0, u_i t_i = 0
    1563: */
    1564: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes,Mat A,Vec f)
    1565: {

    1569:   BCApplyZero_NORTH(dm_stokes,0,A,f);
    1570:   BCApplyZero_NORTH(dm_stokes,1,A,f);
    1571:   BCApplyZero_EAST(dm_stokes,0,A,f);
    1572:   BCApplyZero_SOUTH(dm_stokes,0,A,f);
    1573:   BCApplyZero_SOUTH(dm_stokes,1,A,f);
    1574:   BCApplyZero_WEST(dm_stokes,0,A,f);
    1575:   return(0);
    1576: }

    1578: /*TEST

    1580:    test:
    1581:      suffix: 1
    1582:      args: -no_view
    1583:      requires: !complex double
    1584:      filter: grep -v atomic
    1585:      filter_output: grep -v atomic
    1586:    test:
    1587:      suffix: 1_matis
    1588:      requires: !complex double
    1589:      args: -no_view -dm_mat_type is
    1590:      filter: grep -v atomic
    1591:      filter_output: grep -v atomic
    1592:    testset:
    1593:      nsize: 4
    1594:      requires: !complex double
    1595:      args: -no_view -dm_mat_type is -stokes_ksp_type fetidp -mx 80 -my 80 -stokes_ksp_converged_reason -stokes_ksp_rtol 1.0e-8 -ppcell 2 -nt 4 -randomize_coords -stokes_ksp_error_if_not_converged
    1596:      filter: grep -v atomic
    1597:      filter_output: grep -v atomic
    1598:      test:
    1599:        suffix: fetidp
    1600:        args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
    1601:      test:
    1602:        suffix: fetidp_lumped
    1603:        args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd -stokes_fetidp_pc_lumped -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static
    1604:      test:
    1605:        suffix: fetidp_saddlepoint
    1606:        args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -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
    1607:      test:
    1608:        suffix: fetidp_saddlepoint_lumped
    1609:        args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -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 -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static -stokes_fetidp_pc_lumped
    1610: TEST*/