Actual source code: ex49.c

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

     68: /* Contributed by Dave May */

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

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

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

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

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

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

    103: /*

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

    109:  B = [ d_dx   0   ]
    110:      [  0    d_dy ]
    111:      [ d_dy  d_dx ]

    113:  */

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

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

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

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

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

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

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

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

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

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

    179:   if (det_J) *det_J = J;
    180: }

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

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

    209:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

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

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

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

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

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

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

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

    251:   *_lx = LX;
    252:   *_ly = LY;

    254:   VecDestroy(&vlx);
    255:   VecDestroy(&vly);
    256:   return(0);
    257: }

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

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

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

    292:   PetscFClose(PETSC_COMM_SELF,fp);
    293:   return(0);
    294: }

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

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

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

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

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

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

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

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

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

    357:   PetscFClose(PETSC_COMM_SELF,fp);
    358:   return(0);
    359: }

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

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

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

    390:   DMGetCoordinateDM(da,&cda);
    391:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

    393:   DMCreateLocalVector(da,&local_fields);
    394:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
    395:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
    396:   DMDAVecGetArray(da,local_fields,&_coefficients);

    398:   for (j = sj; j < sj+ny; j++) {
    399:     for (i = si; i < si+nx; i++) {
    400:       PetscScalar coord_x,coord_y;

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

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

    408:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e",
    409:                             PetscRealPart(_coefficients[j][i].E[p]),PetscRealPart(_coefficients[j][i].nu[p]),
    410:                             PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
    411:         PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
    412:       }
    413:     }
    414:   }
    415:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
    416:   VecDestroy(&local_fields);

    418:   PetscFClose(PETSC_COMM_SELF,fp);
    419:   return(0);
    420: }

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

    433:   /* define quadrature rule */
    434:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    436:   /* evaluate integral */
    437:   for (p = 0; p < ngp; p++) {
    438:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
    439:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

    441:     for (i = 0; i < NODES_PER_EL; i++) {
    442:       PetscScalar d_dx_i = GNx_p[0][i];
    443:       PetscScalar d_dy_i = GNx_p[1][i];

    445:       B[0][2*i] = d_dx_i;  B[0][2*i+1] = 0.0;
    446:       B[1][2*i] = 0.0;     B[1][2*i+1] = d_dy_i;
    447:       B[2][2*i] = d_dy_i;  B[2][2*i+1] = d_dx_i;
    448:     }

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

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

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

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

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

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

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

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

    521:   /* node 1 */
    522:   s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0;        /* Ux1 */
    523:   s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1;        /* Uy1 */

    525:   /* node 2 */
    526:   s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0;      /* Ux2 */
    527:   s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1;      /* Uy2 */

    529:   /* node 3 */
    530:   s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0;        /* Ux3 */
    531:   s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1;        /* Uy3 */
    532:   return(0);
    533: }

    535: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
    536: {
    538:   /* get coords for the element */
    539:   el_coords[NSD*0+0] = _coords[ej][ei].x;      el_coords[NSD*0+1] = _coords[ej][ei].y;
    540:   el_coords[NSD*1+0] = _coords[ej+1][ei].x;    el_coords[NSD*1+1] = _coords[ej+1][ei].y;
    541:   el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;  el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
    542:   el_coords[NSD*3+0] = _coords[ej][ei+1].x;    el_coords[NSD*3+1] = _coords[ej][ei+1].y;
    543:   return(0);
    544: }

    546: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
    547: {
    548:   DM                     cda;
    549:   Vec                    coords;
    550:   DMDACoor2d             **_coords;
    551:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    552:   PetscInt               sex,sey,mx,my;
    553:   PetscInt               ei,ej;
    554:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
    555:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    556:   Vec                    local_properties;
    557:   GaussPointCoefficients **props;
    558:   PetscScalar            *prop_E,*prop_nu;
    559:   PetscErrorCode         ierr;

    562:   /* setup for coords */
    563:   DMGetCoordinateDM(elas_da,&cda);
    564:   DMGetCoordinatesLocal(elas_da,&coords);
    565:   DMDAVecGetArray(cda,coords,&_coords);

    567:   /* setup for coefficients */
    568:   DMCreateLocalVector(properties_da,&local_properties);
    569:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
    570:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
    571:   DMDAVecGetArray(properties_da,local_properties,&props);

    573:   DMDAGetElementsCorners(elas_da,&sex,&sey,0);
    574:   DMDAGetElementsSizes(elas_da,&mx,&my,0);
    575:   for (ej = sey; ej < sey+my; ej++) {
    576:     for (ei = sex; ei < sex+mx; ei++) {
    577:       /* get coords for the element */
    578:       GetElementCoords(_coords,ei,ej,el_coords);

    580:       /* get coefficients for the element */
    581:       prop_E  = props[ej][ei].E;
    582:       prop_nu = props[ej][ei].nu;

    584:       /* initialise element stiffness matrix */
    585:       PetscMemzero(Ae,sizeof(Ae));

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

    590:       /* insert element matrix into global matrix */
    591:       DMDAGetElementEqnums_u(u_eqn,ei,ej);
    592:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
    593:     }
    594:   }
    595:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    596:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

    600:   DMDAVecRestoreArray(properties_da,local_properties,&props);
    601:   VecDestroy(&local_properties);
    602:   return(0);
    603: }

    606: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
    607: {
    608:   PetscInt n;

    611:   for (n = 0; n < 4; n++) {
    612:     fields_F[u_eqn[2*n].j][u_eqn[2*n].i].ux_dof     = fields_F[u_eqn[2*n].j][u_eqn[2*n].i].ux_dof+Fe_u[2*n];
    613:     fields_F[u_eqn[2*n+1].j][u_eqn[2*n+1].i].uy_dof = fields_F[u_eqn[2*n+1].j][u_eqn[2*n+1].i].uy_dof+Fe_u[2*n+1];
    614:   }
    615:   return(0);
    616: }

    618: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
    619: {
    620:   DM                     cda;
    621:   Vec                    coords;
    622:   DMDACoor2d             **_coords;
    623:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
    624:   PetscInt               sex,sey,mx,my;
    625:   PetscInt               ei,ej;
    626:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
    627:   PetscScalar            el_coords[NODES_PER_EL*NSD];
    628:   Vec                    local_properties;
    629:   GaussPointCoefficients **props;
    630:   PetscScalar            *prop_fx,*prop_fy;
    631:   Vec                    local_F;
    632:   ElasticityDOF          **ff;
    633:   PetscErrorCode         ierr;

    636:   /* setup for coords */
    637:   DMGetCoordinateDM(elas_da,&cda);
    638:   DMGetCoordinatesLocal(elas_da,&coords);
    639:   DMDAVecGetArray(cda,coords,&_coords);

    641:   /* setup for coefficients */
    642:   DMGetLocalVector(properties_da,&local_properties);
    643:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
    644:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
    645:   DMDAVecGetArray(properties_da,local_properties,&props);

    647:   /* get acces to the vector */
    648:   DMGetLocalVector(elas_da,&local_F);
    649:   VecZeroEntries(local_F);
    650:   DMDAVecGetArray(elas_da,local_F,&ff);

    652:   DMDAGetElementsCorners(elas_da,&sex,&sey,0);
    653:   DMDAGetElementsSizes(elas_da,&mx,&my,0);
    654:   for (ej = sey; ej < sey+my; ej++) {
    655:     for (ei = sex; ei < sex+mx; ei++) {
    656:       /* get coords for the element */
    657:       GetElementCoords(_coords,ei,ej,el_coords);

    659:       /* get coefficients for the element */
    660:       prop_fx = props[ej][ei].fx;
    661:       prop_fy = props[ej][ei].fy;

    663:       /* initialise element stiffness matrix */
    664:       PetscMemzero(Fe,sizeof(Fe));

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

    669:       /* insert element matrix into global matrix */
    670:       DMDAGetElementEqnums_u(u_eqn,ei,ej);

    672:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
    673:     }
    674:   }

    676:   DMDAVecRestoreArray(elas_da,local_F,&ff);
    677:   DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
    678:   DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
    679:   DMRestoreLocalVector(elas_da,&local_F);

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

    683:   DMDAVecRestoreArray(properties_da,local_properties,&props);
    684:   DMRestoreLocalVector(properties_da,&local_properties);
    685:   return(0);
    686: }

    688: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
    689: {
    690:   DM                     elas_da,da_prop;
    691:   PetscInt               u_dof,dof,stencil_width;
    692:   Mat                    A;
    693:   PetscInt               mxl,myl;
    694:   DM                     prop_cda,vel_cda;
    695:   Vec                    prop_coords,vel_coords;
    696:   PetscInt               si,sj,nx,ny,i,j,p;
    697:   Vec                    f,X;
    698:   PetscInt               prop_dof,prop_stencil_width;
    699:   Vec                    properties,l_properties;
    700:   MatNullSpace           matnull;
    701:   PetscReal              dx,dy;
    702:   PetscInt               M,N;
    703:   DMDACoor2d             **_prop_coords,**_vel_coords;
    704:   GaussPointCoefficients **element_props;
    705:   KSP                    ksp_E;
    706:   PetscInt               coefficient_structure = 0;
    707:   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
    708:   PetscBool              use_gp_coords = PETSC_FALSE;
    709:   PetscBool              use_nonsymbc  = PETSC_FALSE;
    710:   PetscBool              no_view       = PETSC_FALSE;
    711:   PetscBool              flg;
    712:   PetscErrorCode         ierr;

    715:   /* Generate the da for velocity and pressure */
    716:   /*
    717:    We use Q1 elements for the temperature.
    718:    FEM has a 9-point stencil (BOX) or connectivity pattern
    719:    Num nodes in each direction is mx+1, my+1
    720:    */
    721:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
    722:   dof           = u_dof;
    723:   stencil_width = 1;
    724:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&elas_da);

    726:   DMSetMatType(elas_da,MATAIJ);
    727:   DMSetFromOptions(elas_da);
    728:   DMSetUp(elas_da);

    730:   DMDASetFieldName(elas_da,0,"Ux");
    731:   DMDASetFieldName(elas_da,1,"Uy");

    733:   /* unit box [0,1] x [0,1] */
    734:   DMDASetUniformCoordinates(elas_da,0.0,1.0,0.0,1.0,0.0,1.0);

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

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

    744:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
    745:   prop_stencil_width = 0;
    746:   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);
    747:   DMSetFromOptions(da_prop);
    748:   DMSetUp(da_prop);

    750:   PetscFree(lx);
    751:   PetscFree(ly);

    753:   /* define centroid positions */
    754:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
    755:   dx   = 1.0/((PetscReal)(M));
    756:   dy   = 1.0/((PetscReal)(N));

    758:   DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,0.0,1.0);

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

    763:   DMCreateGlobalVector(da_prop,&properties);
    764:   DMCreateLocalVector(da_prop,&l_properties);
    765:   DMDAVecGetArray(da_prop,l_properties,&element_props);

    767:   DMGetCoordinateDM(da_prop,&prop_cda);
    768:   DMGetCoordinatesLocal(da_prop,&prop_coords);
    769:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

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

    773:   DMGetCoordinateDM(elas_da,&vel_cda);
    774:   DMGetCoordinatesLocal(elas_da,&vel_coords);
    775:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);

    777:   /* interpolate the coordinates */
    778:   for (j = sj; j < sj+ny; j++) {
    779:     for (i = si; i < si+nx; i++) {
    780:       PetscInt    ngp;
    781:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
    782:       PetscScalar el_coords[8];

    784:       GetElementCoords(_vel_coords,i,j,el_coords);
    785:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

    787:       for (p = 0; p < GAUSS_POINTS; p++) {
    788:         PetscScalar gp_x,gp_y;
    789:         PetscInt    n;
    790:         PetscScalar xi_p[2],Ni_p[4];

    792:         xi_p[0] = gp_xi[p][0];
    793:         xi_p[1] = gp_xi[p][1];
    794:         ConstructQ12D_Ni(xi_p,Ni_p);

    796:         gp_x = 0.0;
    797:         gp_y = 0.0;
    798:         for (n = 0; n < NODES_PER_EL; n++) {
    799:           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
    800:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
    801:         }
    802:         element_props[j][i].gp_coords[2*p]   = gp_x;
    803:         element_props[j][i].gp_coords[2*p+1] = gp_y;
    804:       }
    805:     }
    806:   }

    808:   /* define the coefficients */
    809:   PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,&flg);

    811:   for (j = sj; j < sj+ny; j++) {
    812:     for (i = si; i < si+nx; i++) {
    813:       PetscScalar              centroid_x = _prop_coords[j][i].x; /* centroids of cell */
    814:       PetscScalar              centroid_y = _prop_coords[j][i].y;
    815:       PETSC_UNUSED PetscScalar coord_x,coord_y;

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

    820:         opts_E  = 1.0;
    821:         opts_nu = 0.33;
    822:         PetscOptionsGetScalar(NULL,NULL,"-iso_E",&opts_E,&flg);
    823:         PetscOptionsGetScalar(NULL,NULL,"-iso_nu",&opts_nu,&flg);

    825:         for (p = 0; p < GAUSS_POINTS; p++) {
    826:           element_props[j][i].E[p]  = opts_E;
    827:           element_props[j][i].nu[p] = opts_nu;

    829:           element_props[j][i].fx[p] = 0.0;
    830:           element_props[j][i].fy[p] = 0.0;
    831:         }
    832:       } else if (coefficient_structure == 1) { /* step */
    833:         PetscScalar opts_E0,opts_nu0,opts_xc;
    834:         PetscScalar opts_E1,opts_nu1;

    836:         opts_E0  = opts_E1  = 1.0;
    837:         opts_nu0 = opts_nu1 = 0.333;
    838:         opts_xc  = 0.5;
    839:         PetscOptionsGetScalar(NULL,NULL,"-step_E0",&opts_E0,&flg);
    840:         PetscOptionsGetScalar(NULL,NULL,"-step_nu0",&opts_nu0,&flg);
    841:         PetscOptionsGetScalar(NULL,NULL,"-step_E1",&opts_E1,&flg);
    842:         PetscOptionsGetScalar(NULL,NULL,"-step_nu1",&opts_nu1,&flg);
    843:         PetscOptionsGetScalar(NULL,NULL,"-step_xc",&opts_xc,&flg);

    845:         for (p = 0; p < GAUSS_POINTS; p++) {
    846:           coord_x = centroid_x;
    847:           coord_y = centroid_y;
    848:           if (use_gp_coords) {
    849:             coord_x = element_props[j][i].gp_coords[2*p];
    850:             coord_y = element_props[j][i].gp_coords[2*p+1];
    851:           }

    853:           element_props[j][i].E[p]  = opts_E0;
    854:           element_props[j][i].nu[p] = opts_nu0;
    855:           if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
    856:             element_props[j][i].E[p]  = opts_E1;
    857:             element_props[j][i].nu[p] = opts_nu1;
    858:           }

    860:           element_props[j][i].fx[p] = 0.0;
    861:           element_props[j][i].fy[p] = 0.0;
    862:         }
    863:       } else if (coefficient_structure == 2) { /* brick */
    864:         PetscReal values_E[10];
    865:         PetscReal values_nu[10];
    866:         PetscInt  nbricks,maxnbricks;
    867:         PetscInt  index,span;
    868:         PetscInt  jj;

    870:         flg        = PETSC_FALSE;
    871:         maxnbricks = 10;
    872:         PetscOptionsGetRealArray(NULL,NULL, "-brick_E",values_E,&maxnbricks,&flg);
    873:         nbricks    = maxnbricks;
    874:         if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of E values for each brick");

    876:         flg        = PETSC_FALSE;
    877:         maxnbricks = 10;
    878:         PetscOptionsGetRealArray(NULL,NULL, "-brick_nu",values_nu,&maxnbricks,&flg);
    879:         if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of nu values for each brick");
    880:         if (maxnbricks != nbricks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply equal numbers of values for E and nu");

    882:         span = 1;
    883:         PetscOptionsGetInt(NULL,NULL,"-brick_span",&span,&flg);

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

    890:         for (p = 0; p < GAUSS_POINTS; p++) {
    891:           element_props[j][i].E[p]  = values_E[index];
    892:           element_props[j][i].nu[p] = values_nu[index];
    893:         }
    894:       } else if (coefficient_structure == 3) { /* sponge */
    895:         PetscScalar opts_E0,opts_nu0;
    896:         PetscScalar opts_E1,opts_nu1;
    897:         PetscInt    opts_t,opts_w;
    898:         PetscInt    ii,jj,ci,cj;

    900:         opts_E0  = opts_E1  = 1.0;
    901:         opts_nu0 = opts_nu1 = 0.333;
    902:         PetscOptionsGetScalar(NULL,NULL,"-sponge_E0",&opts_E0,&flg);
    903:         PetscOptionsGetScalar(NULL,NULL,"-sponge_nu0",&opts_nu0,&flg);
    904:         PetscOptionsGetScalar(NULL,NULL,"-sponge_E1",&opts_E1,&flg);
    905:         PetscOptionsGetScalar(NULL,NULL,"-sponge_nu1",&opts_nu1,&flg);

    907:         opts_t = opts_w = 1;
    908:         PetscOptionsGetInt(NULL,NULL,"-sponge_t",&opts_t,&flg);
    909:         PetscOptionsGetInt(NULL,NULL,"-sponge_w",&opts_w,&flg);

    911:         ii = (i)/(opts_t+opts_w+opts_t);
    912:         jj = (j)/(opts_t+opts_w+opts_t);

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

    917:         for (p = 0; p < GAUSS_POINTS; p++) {
    918:           element_props[j][i].E[p]  = opts_E0;
    919:           element_props[j][i].nu[p] = opts_nu0;
    920:         }
    921:         if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
    922:           if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
    923:             for (p = 0; p < GAUSS_POINTS; p++) {
    924:               element_props[j][i].E[p]  = opts_E1;
    925:               element_props[j][i].nu[p] = opts_nu1;
    926:             }
    927:           }
    928:         }
    929:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
    930:     }
    931:   }
    932:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

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

    936:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
    937:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
    938:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);

    940:   PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
    941:   if (!no_view) {
    942:     DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
    943:     DMDACoordViewGnuplot2d(elas_da,"mesh");
    944:   }

    946:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
    947:   DMCreateMatrix(elas_da,&A);
    948:   DMGetCoordinates(elas_da,&vel_coords);
    949:   MatNullSpaceCreateRigidBody(vel_coords,&matnull);
    950:   MatSetNearNullSpace(A,matnull);
    951:   MatNullSpaceDestroy(&matnull);
    952:   MatCreateVecs(A,&f,&X);

    954:   /* assemble A11 */
    955:   MatZeroEntries(A);
    956:   VecZeroEntries(f);

    958:   AssembleA_Elasticity(A,elas_da,da_prop,properties);
    959:   /* build force vector */
    960:   AssembleF_Elasticity(f,elas_da,da_prop,properties);

    962:   KSPCreate(PETSC_COMM_WORLD,&ksp_E);
    963:   KSPSetOptionsPrefix(ksp_E,"elas_");  /* elasticity */

    965:   PetscOptionsGetBool(NULL,NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
    966:   /* solve */
    967:   if (!use_nonsymbc) {
    968:     Mat        AA;
    969:     Vec        ff,XX;
    970:     IS         is;
    971:     VecScatter scat;

    973:     DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
    974:     VecDuplicate(ff,&XX);

    976:     KSPSetOperators(ksp_E,AA,AA);
    977:     KSPSetFromOptions(ksp_E);

    979:     KSPSolve(ksp_E,ff,XX);

    981:     /* push XX back into X */
    982:     DMDABCApplyCompression(elas_da,NULL,X);

    984:     VecScatterCreate(XX,NULL,X,is,&scat);
    985:     VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
    986:     VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
    987:     VecScatterDestroy(&scat);

    989:     MatDestroy(&AA);
    990:     VecDestroy(&ff);
    991:     VecDestroy(&XX);
    992:     ISDestroy(&is);
    993:   } else {
    994:     DMDABCApplyCompression(elas_da,A,f);

    996:     KSPSetOperators(ksp_E,A,A);
    997:     KSPSetFromOptions(ksp_E);

    999:     KSPSolve(ksp_E,f,X);
    1000:   }

    1002:   if (!no_view) {DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");}
    1003:   KSPDestroy(&ksp_E);

    1005:   VecDestroy(&X);
    1006:   VecDestroy(&f);
    1007:   MatDestroy(&A);

    1009:   DMDestroy(&elas_da);
    1010:   DMDestroy(&da_prop);

    1012:   VecDestroy(&properties);
    1013:   VecDestroy(&l_properties);
    1014:   return(0);
    1015: }

    1017: int main(int argc,char **args)
    1018: {
    1020:   PetscInt       mx,my;

    1022:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
    1023:   mx   = my = 10;
    1024:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
    1025:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
    1026:   solve_elasticity_2d(mx,my);
    1027:   PetscFinalize();
    1028:   return ierr;
    1029: }

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

    1033: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
    1034: {
    1035:   DM                     cda;
    1036:   Vec                    coords;
    1037:   PetscInt               si,sj,nx,ny,i,j;
    1038:   PetscInt               M,N;
    1039:   DMDACoor2d             **_coords;
    1040:   const PetscInt         *g_idx;
    1041:   PetscInt               *bc_global_ids;
    1042:   PetscScalar            *bc_vals;
    1043:   PetscInt               nbcs;
    1044:   PetscInt               n_dofs;
    1045:   PetscErrorCode         ierr;
    1046:   ISLocalToGlobalMapping ltogm;

    1049:   /* enforce bc's */
    1050:   DMGetLocalToGlobalMapping(da,&ltogm);
    1051:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1053:   DMGetCoordinateDM(da,&cda);
    1054:   DMGetCoordinatesLocal(da,&coords);
    1055:   DMDAVecGetArray(cda,coords,&_coords);
    1056:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1057:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1059:   /* --- */

    1061:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
    1062:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

    1067:   i = nx-1;
    1068:   for (j = 0; j < ny; j++) {
    1069:     PetscInt                 local_id;
    1070:     PETSC_UNUSED PetscScalar coordx,coordy;

    1072:     local_id = i+j*nx;

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

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

    1079:     bc_vals[j] =  bc_val;
    1080:   }
    1081:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1082:   nbcs = 0;
    1083:   if ((si+nx) == (M)) nbcs = ny;

    1085:   if (b) {
    1086:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1087:     VecAssemblyBegin(b);
    1088:     VecAssemblyEnd(b);
    1089:   }
    1090:   if (A) {
    1091:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
    1092:   }

    1094:   PetscFree(bc_vals);
    1095:   PetscFree(bc_global_ids);

    1097:   DMDAVecRestoreArray(cda,coords,&_coords);
    1098:   return(0);
    1099: }

    1101: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
    1102: {
    1103:   DM                     cda;
    1104:   Vec                    coords;
    1105:   PetscInt               si,sj,nx,ny,i,j;
    1106:   PetscInt               M,N;
    1107:   DMDACoor2d             **_coords;
    1108:   const PetscInt         *g_idx;
    1109:   PetscInt               *bc_global_ids;
    1110:   PetscScalar            *bc_vals;
    1111:   PetscInt               nbcs;
    1112:   PetscInt               n_dofs;
    1113:   PetscErrorCode         ierr;
    1114:   ISLocalToGlobalMapping ltogm;

    1117:   /* enforce bc's */
    1118:   DMGetLocalToGlobalMapping(da,&ltogm);
    1119:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

    1121:   DMGetCoordinateDM(da,&cda);
    1122:   DMGetCoordinatesLocal(da,&coords);
    1123:   DMDAVecGetArray(cda,coords,&_coords);
    1124:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
    1125:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

    1127:   /* --- */

    1129:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
    1130:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

    1135:   i = 0;
    1136:   for (j = 0; j < ny; j++) {
    1137:     PetscInt                 local_id;
    1138:     PETSC_UNUSED PetscScalar coordx,coordy;

    1140:     local_id = i+j*nx;

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

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

    1147:     bc_vals[j] =  bc_val;
    1148:   }
    1149:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
    1150:   nbcs = 0;
    1151:   if (si == 0) nbcs = ny;

    1153:   if (b) {
    1154:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
    1155:     VecAssemblyBegin(b);
    1156:     VecAssemblyEnd(b);
    1157:   }
    1158:   if (A) {
    1159:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
    1160:   }

    1162:   PetscFree(bc_vals);
    1163:   PetscFree(bc_global_ids);

    1165:   DMDAVecRestoreArray(cda,coords,&_coords);
    1166:   return(0);
    1167: }

    1169: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
    1170: {

    1174:   BCApply_EAST(elas_da,0,-1.0,A,f);
    1175:   BCApply_EAST(elas_da,1, 0.0,A,f);
    1176:   BCApply_WEST(elas_da,0,1.0,A,f);
    1177:   BCApply_WEST(elas_da,1,0.0,A,f);
    1178:   return(0);
    1179: }

    1181: static PetscErrorCode Orthogonalize(PetscInt n,Vec *vecs)
    1182: {
    1183:   PetscInt       i,j;
    1184:   PetscScalar    dot;

    1188:   for (i=0; i<n; i++) {
    1189:   VecNormalize(vecs[i],NULL);
    1190:      for (j=i+1; j<n; j++) {
    1191:        VecDot(vecs[i],vecs[j],&dot);
    1192:        VecAXPY(vecs[j],-dot,vecs[i]);
    1193:      }
    1194:   }
    1195:   return(0);
    1196: }

    1198: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
    1199: {
    1201:   PetscInt       start,end,m;
    1202:   PetscInt       *unconstrained;
    1203:   PetscInt       cnt,i;
    1204:   Vec            x;
    1205:   PetscScalar    *_x;
    1206:   IS             is;
    1207:   VecScatter     scat;

    1210:   /* push bc's into f and A */
    1211:   VecDuplicate(f,&x);
    1212:   BCApply_EAST(elas_da,0,-1.0,A,x);
    1213:   BCApply_EAST(elas_da,1, 0.0,A,x);
    1214:   BCApply_WEST(elas_da,0,1.0,A,x);
    1215:   BCApply_WEST(elas_da,1,0.0,A,x);

    1217:   /* define which dofs are not constrained */
    1218:   VecGetLocalSize(x,&m);
    1219:   PetscMalloc1(m,&unconstrained);
    1220:   VecGetOwnershipRange(x,&start,&end);
    1221:   VecGetArray(x,&_x);
    1222:   cnt  = 0;
    1223:   for (i = 0; i < m; i+=2) {
    1224:     PetscReal val1,val2;

    1226:     val1 = PetscRealPart(_x[i]);
    1227:     val2 = PetscRealPart(_x[i+1]);
    1228:     if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
    1229:       unconstrained[cnt] = start + i;
    1230:       cnt++;
    1231:       unconstrained[cnt] = start + i + 1;
    1232:       cnt++;
    1233:     }
    1234:   }
    1235:   VecRestoreArray(x,&_x);

    1237:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
    1238:   PetscFree(unconstrained);
    1239:   ISSetBlockSize(is,2);

    1241:   /* define correction for dirichlet in the rhs */
    1242:   MatMult(A,x,f);
    1243:   VecScale(f,-1.0);

    1245:   /* get new matrix */
    1246:   MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
    1247:   /* get new vector */
    1248:   MatCreateVecs(*AA,NULL,ff);

    1250:   VecScatterCreate(f,is,*ff,NULL,&scat);
    1251:   VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
    1252:   VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);

    1254:   {                             /* Constrain near-null space */
    1255:     PetscInt     nvecs;
    1256:     const        Vec *vecs;
    1257:     Vec          *uvecs;
    1258:     PetscBool    has_const;
    1259:     MatNullSpace mnull,unull;

    1261:     MatGetNearNullSpace(A,&mnull);
    1262:     MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);
    1263:     VecDuplicateVecs(*ff,nvecs,&uvecs);
    1264:     for (i=0; i<nvecs; i++) {
    1265:       VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
    1266:       VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
    1267:     }
    1268:     Orthogonalize(nvecs,uvecs);
    1269:     MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,nvecs,uvecs,&unull);
    1270:     MatSetNearNullSpace(*AA,unull);
    1271:     MatNullSpaceDestroy(&unull);
    1272:     VecDestroyVecs(nvecs,&uvecs);
    1273:   }

    1275:   VecScatterDestroy(&scat);

    1277:   *dofs = is;
    1278:   VecDestroy(&x);
    1279:   return(0);
    1280: }

    1283: /*TEST

    1285:    build:
    1286:       requires: !complex !single

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

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

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

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

    1308:    test:
    1309:       suffix: 5
    1310:       nsize: 3
    1311:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type gamg -elas_mg_levels_ksp_type chebyshev -elas_mg_levels_ksp_max_it 1 -elas_mg_levels_ksp_chebyshev_esteig 0.2,1.1 -elas_mg_levels_pc_type jacobi

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

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

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

    1328:    test:
    1329:       suffix: hypre_nullspace
    1330:       requires: hypre
    1331:       args: -elas_ksp_monitor_short -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -c_str 2 -brick_E 1,1e-6,1e-2 -brick_nu .3,.2,.4 -brick_span 8 -elas_pc_type hypre -elas_pc_hypre_boomeramg_nodal_coarsen 6 -elas_pc_hypre_boomeramg_vec_interp_variant 3 -elas_ksp_view

    1333:    test:
    1334:       nsize: 4
    1335:       suffix: bddc
    1336:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_pc_type bddc -elas_pc_bddc_monolithic

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

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

    1348:    test:
    1349:       nsize: 4
    1350:       suffix: fetidp_unsym_deluxe
    1351:       args: -elas_ksp_monitor_short -no_view -elas_ksp_converged_reason -elas_ksp_type fetidp -elas_fetidp_ksp_type cg -elas_ksp_norm_type natural -mx 22 -my 22 -dm_mat_type is -elas_fetidp_bddc_pc_bddc_monolithic -use_nonsymbc -elas_fetidp_bddc_pc_bddc_use_deluxe_scaling -elas_fetidp_bddc_sub_schurs_symmetric 0 -elas_fetidp_bddc_pc_bddc_deluxe_singlemat

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

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

    1363:    test:
    1364:       nsize: 9
    1365:       suffix: bddc_subdomainjump_deluxe
    1366:       args: -mx 30 -my 30 -elas_ksp_monitor_short -no_view -c_str 2 -brick_E 1,1000 -brick_nu 0.4,0.2 -brick_span 10  -elas_ksp_type cg -elas_pc_type bddc -elas_pc_bddc_monolithic -dm_mat_type is -elas_pc_bddc_use_deluxe_scaling -elas_ksp_norm_type natural -elas_pc_bddc_schur_layers 1
    1367: TEST*/