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

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

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

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

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

101: /*

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

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

111:  */

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

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

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

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

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

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

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

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

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

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

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

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

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

205:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

209:   proc_J = rank/cpu_x;
210:   proc_I = rank-cpu_x*proc_J;

212:   PetscMalloc1(cpu_x,&LX);
213:   PetscMalloc1(cpu_y,&LY);

215:   DMDAGetElementsSizes(da,&local_mx,&local_my,NULL);
216:   VecCreate(PETSC_COMM_WORLD,&vlx);
217:   VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
218:   VecSetFromOptions(vlx);

220:   VecCreate(PETSC_COMM_WORLD,&vly);
221:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
222:   VecSetFromOptions(vly);

224:   VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
225:   VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
226:   VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
227:   VecAssemblyBegin(vly);VecAssemblyEnd(vly);

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

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

247:   *_lx = LX;
248:   *_ly = LY;

250:   VecDestroy(&vlx);
251:   VecDestroy(&vly);
252:   return 0;
253: }

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

266:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
267:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
268:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
270:   PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);

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

287:   PetscFClose(PETSC_COMM_SELF,fp);
288:   return 0;
289: }

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

305:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
306:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
307:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);

310:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
311:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
312:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
313:   for (d = 0; d < n_dofs; d++) {
314:     DMDAGetFieldName(da,d,&field_name);
315:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
316:   }
317:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");

319:   DMGetCoordinateDM(da,&cda);
320:   DMGetCoordinatesLocal(da,&coords);
321:   DMDAVecGetArray(cda,coords,&_coords);
322:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

324:   DMCreateLocalVector(da,&local_fields);
325:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
326:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
327:   VecGetArray(local_fields,&_fields);

329:   for (j = sj; j < sj+ny; j++) {
330:     for (i = si; i < si+nx; i++) {
331:       PetscScalar coord_x,coord_y;
332:       PetscScalar field_d;

334:       coord_x = _coords[j][i].x;
335:       coord_y = _coords[j][i].y;

337:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
338:       for (d = 0; d < n_dofs; d++) {
339:         field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
340:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
341:       }
342:       PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
343:     }
344:   }
345:   VecRestoreArray(local_fields,&_fields);
346:   VecDestroy(&local_fields);

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

350:   PetscFClose(PETSC_COMM_SELF,fp);
351:   return 0;
352: }

354: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
355: {
356:   DM                     cda;
357:   Vec                    local_fields;
358:   FILE                   *fp;
359:   char                   fname[PETSC_MAX_PATH_LEN];
360:   const char             *field_name;
361:   PetscMPIInt            rank;
362:   PetscInt               si,sj,nx,ny,i,j,p;
363:   PetscInt               n_dofs,d;
364:   GaussPointCoefficients **_coefficients;

367:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
368:   PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
369:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);

372:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
373:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
374:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
375:   for (d = 0; d < n_dofs; d++) {
376:     DMDAGetFieldName(da,d,&field_name);
377:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
378:   }
379:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");

381:   DMGetCoordinateDM(da,&cda);
382:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

384:   DMCreateLocalVector(da,&local_fields);
385:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
386:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
387:   DMDAVecGetArray(da,local_fields,&_coefficients);

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

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

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

399:         PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e\n",
400:                                (double)PetscRealPart(_coefficients[j][i].E[p]),(double)PetscRealPart(_coefficients[j][i].nu[p]),
401:                                (double)PetscRealPart(_coefficients[j][i].fx[p]),(double)PetscRealPart(_coefficients[j][i].fy[p])));
402:       }
403:     }
404:   }
405:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
406:   VecDestroy(&local_fields);

408:   PetscFClose(PETSC_COMM_SELF,fp);
409:   return 0;
410: }

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

423:   /* define quadrature rule */
424:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

426:   /* evaluate integral */
427:   for (p = 0; p < ngp; p++) {
428:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
429:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

431:     for (i = 0; i < NODES_PER_EL; i++) {
432:       PetscScalar d_dx_i = GNx_p[0][i];
433:       PetscScalar d_dy_i = GNx_p[1][i];

435:       B[0][2*i] = d_dx_i;  B[0][2*i+1] = 0.0;
436:       B[1][2*i] = 0.0;     B[1][2*i+1] = d_dy_i;
437:       B[2][2*i] = d_dy_i;  B[2][2*i+1] = d_dx_i;
438:     }

440:     /* form D for the quadrature point */
441:     prop_E          = E[p];
442:     prop_nu         = nu[p];
443:     factor          = prop_E / ((1.0+prop_nu)*(1.0-2.0*prop_nu));
444:     constit_D[0][0] = 1.0-prop_nu;  constit_D[0][1] = prop_nu;      constit_D[0][2] = 0.0;
445:     constit_D[1][0] = prop_nu;      constit_D[1][1] = 1.0-prop_nu;  constit_D[1][2] = 0.0;
446:     constit_D[2][0] = 0.0;          constit_D[2][1] = 0.0;          constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
447:     for (i = 0; i < 3; i++) {
448:       for (j = 0; j < 3; j++) {
449:         constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
450:       }
451:     }

453:     /* form Bt tildeD B */
454:     /*
455:      Ke_ij = Bt_ik . D_kl . B_lj
456:      = B_ki . D_kl . B_lj
457:      */
458:     for (i = 0; i < 8; i++) {
459:       for (j = 0; j < 8; j++) {
460:         for (k = 0; k < 3; k++) {
461:           for (l = 0; l < 3; l++) {
462:             Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
463:           }
464:         }
465:       }
466:     }

468:   } /* end quadrature */
469: }

471: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
472: {
473:   PetscInt    ngp;
474:   PetscScalar gp_xi[GAUSS_POINTS][2];
475:   PetscScalar gp_weight[GAUSS_POINTS];
476:   PetscInt    p,i;
477:   PetscScalar Ni_p[NODES_PER_EL];
478:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
479:   PetscScalar J_p,fac;

481:   /* define quadrature rule */
482:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

484:   /* evaluate integral */
485:   for (p = 0; p < ngp; p++) {
486:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
487:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
488:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
489:     fac = gp_weight[p]*J_p;

491:     for (i = 0; i < NODES_PER_EL; i++) {
492:       Fe[NSD*i]   += fac*Ni_p[i]*fx[p];
493:       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
494:     }
495:   }
496: }

498: /*
499:  i,j are the element indices
500:  The unknown is a vector quantity.
501:  The s[].c is used to indicate the degree of freedom.
502:  */
503: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)
504: {
506:   /* displacement */
507:   /* node 0 */
508:   s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0;          /* Ux0 */
509:   s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1;          /* Uy0 */

511:   /* node 1 */
512:   s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0;        /* Ux1 */
513:   s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1;        /* Uy1 */

515:   /* node 2 */
516:   s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0;      /* Ux2 */
517:   s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1;      /* Uy2 */

519:   /* node 3 */
520:   s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0;        /* Ux3 */
521:   s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1;        /* Uy3 */
522:   return 0;
523: }

525: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
526: {
528:   /* get coords for the element */
529:   el_coords[NSD*0+0] = _coords[ej][ei].x;      el_coords[NSD*0+1] = _coords[ej][ei].y;
530:   el_coords[NSD*1+0] = _coords[ej+1][ei].x;    el_coords[NSD*1+1] = _coords[ej+1][ei].y;
531:   el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;  el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
532:   el_coords[NSD*3+0] = _coords[ej][ei+1].x;    el_coords[NSD*3+1] = _coords[ej][ei+1].y;
533:   return 0;
534: }

536: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
537: {
538:   DM                     cda;
539:   Vec                    coords;
540:   DMDACoor2d             **_coords;
541:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
542:   PetscInt               sex,sey,mx,my;
543:   PetscInt               ei,ej;
544:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
545:   PetscScalar            el_coords[NODES_PER_EL*NSD];
546:   Vec                    local_properties;
547:   GaussPointCoefficients **props;
548:   PetscScalar            *prop_E,*prop_nu;

551:   /* setup for coords */
552:   DMGetCoordinateDM(elas_da,&cda);
553:   DMGetCoordinatesLocal(elas_da,&coords);
554:   DMDAVecGetArray(cda,coords,&_coords);

556:   /* setup for coefficients */
557:   DMCreateLocalVector(properties_da,&local_properties);
558:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
559:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
560:   DMDAVecGetArray(properties_da,local_properties,&props);

562:   DMDAGetElementsCorners(elas_da,&sex,&sey,0);
563:   DMDAGetElementsSizes(elas_da,&mx,&my,0);
564:   for (ej = sey; ej < sey+my; ej++) {
565:     for (ei = sex; ei < sex+mx; ei++) {
566:       /* get coords for the element */
567:       GetElementCoords(_coords,ei,ej,el_coords);

569:       /* get coefficients for the element */
570:       prop_E  = props[ej][ei].E;
571:       prop_nu = props[ej][ei].nu;

573:       /* initialise element stiffness matrix */
574:       PetscMemzero(Ae,sizeof(Ae));

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

579:       /* insert element matrix into global matrix */
580:       DMDAGetElementEqnums_u(u_eqn,ei,ej);
581:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
582:     }
583:   }
584:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
585:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

589:   DMDAVecRestoreArray(properties_da,local_properties,&props);
590:   VecDestroy(&local_properties);
591:   return 0;
592: }

594: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
595: {
596:   PetscInt n;

599:   for (n = 0; n < 4; n++) {
600:     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];
601:     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];
602:   }
603:   return 0;
604: }

606: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
607: {
608:   DM                     cda;
609:   Vec                    coords;
610:   DMDACoor2d             **_coords;
611:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
612:   PetscInt               sex,sey,mx,my;
613:   PetscInt               ei,ej;
614:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
615:   PetscScalar            el_coords[NODES_PER_EL*NSD];
616:   Vec                    local_properties;
617:   GaussPointCoefficients **props;
618:   PetscScalar            *prop_fx,*prop_fy;
619:   Vec                    local_F;
620:   ElasticityDOF          **ff;

623:   /* setup for coords */
624:   DMGetCoordinateDM(elas_da,&cda);
625:   DMGetCoordinatesLocal(elas_da,&coords);
626:   DMDAVecGetArray(cda,coords,&_coords);

628:   /* setup for coefficients */
629:   DMGetLocalVector(properties_da,&local_properties);
630:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
631:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
632:   DMDAVecGetArray(properties_da,local_properties,&props);

634:   /* get access to the vector */
635:   DMGetLocalVector(elas_da,&local_F);
636:   VecZeroEntries(local_F);
637:   DMDAVecGetArray(elas_da,local_F,&ff);

639:   DMDAGetElementsCorners(elas_da,&sex,&sey,0);
640:   DMDAGetElementsSizes(elas_da,&mx,&my,0);
641:   for (ej = sey; ej < sey+my; ej++) {
642:     for (ei = sex; ei < sex+mx; ei++) {
643:       /* get coords for the element */
644:       GetElementCoords(_coords,ei,ej,el_coords);

646:       /* get coefficients for the element */
647:       prop_fx = props[ej][ei].fx;
648:       prop_fy = props[ej][ei].fy;

650:       /* initialise element stiffness matrix */
651:       PetscMemzero(Fe,sizeof(Fe));

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

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

659:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
660:     }
661:   }

663:   DMDAVecRestoreArray(elas_da,local_F,&ff);
664:   DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
665:   DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
666:   DMRestoreLocalVector(elas_da,&local_F);

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

670:   DMDAVecRestoreArray(properties_da,local_properties,&props);
671:   DMRestoreLocalVector(properties_da,&local_properties);
672:   return 0;
673: }

675: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
676: {
677:   DM                     elas_da,da_prop;
678:   PetscInt               u_dof,dof,stencil_width;
679:   Mat                    A;
680:   PetscInt               mxl,myl;
681:   DM                     prop_cda,vel_cda;
682:   Vec                    prop_coords,vel_coords;
683:   PetscInt               si,sj,nx,ny,i,j,p;
684:   Vec                    f,X;
685:   PetscInt               prop_dof,prop_stencil_width;
686:   Vec                    properties,l_properties;
687:   MatNullSpace           matnull;
688:   PetscReal              dx,dy;
689:   PetscInt               M,N;
690:   DMDACoor2d             **_prop_coords,**_vel_coords;
691:   GaussPointCoefficients **element_props;
692:   KSP                    ksp_E;
693:   PetscInt               coefficient_structure = 0;
694:   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
695:   PetscBool              use_gp_coords = PETSC_FALSE;
696:   PetscBool              use_nonsymbc  = PETSC_FALSE;
697:   PetscBool              no_view       = PETSC_FALSE;
698:   PetscBool              flg;

701:   /* Generate the da for velocity and pressure */
702:   /*
703:    We use Q1 elements for the temperature.
704:    FEM has a 9-point stencil (BOX) or connectivity pattern
705:    Num nodes in each direction is mx+1, my+1
706:    */
707:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
708:   dof           = u_dof;
709:   stencil_width = 1;
710:   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);

712:   DMSetMatType(elas_da,MATAIJ);
713:   DMSetFromOptions(elas_da);
714:   DMSetUp(elas_da);

716:   DMDASetFieldName(elas_da,0,"Ux");
717:   DMDASetFieldName(elas_da,1,"Uy");

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

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

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

730:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
731:   prop_stencil_width = 0;
732:   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);
733:   DMSetFromOptions(da_prop);
734:   DMSetUp(da_prop);

736:   PetscFree(lx);
737:   PetscFree(ly);

739:   /* define centroid positions */
740:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
741:   dx   = 1.0/((PetscReal)(M));
742:   dy   = 1.0/((PetscReal)(N));

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

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

749:   DMCreateGlobalVector(da_prop,&properties);
750:   DMCreateLocalVector(da_prop,&l_properties);
751:   DMDAVecGetArray(da_prop,l_properties,&element_props);

753:   DMGetCoordinateDM(da_prop,&prop_cda);
754:   DMGetCoordinatesLocal(da_prop,&prop_coords);
755:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

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

759:   DMGetCoordinateDM(elas_da,&vel_cda);
760:   DMGetCoordinatesLocal(elas_da,&vel_coords);
761:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);

763:   /* interpolate the coordinates */
764:   for (j = sj; j < sj+ny; j++) {
765:     for (i = si; i < si+nx; i++) {
766:       PetscInt    ngp;
767:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
768:       PetscScalar el_coords[8];

770:       GetElementCoords(_vel_coords,i,j,el_coords);
771:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

773:       for (p = 0; p < GAUSS_POINTS; p++) {
774:         PetscScalar gp_x,gp_y;
775:         PetscInt    n;
776:         PetscScalar xi_p[2],Ni_p[4];

778:         xi_p[0] = gp_xi[p][0];
779:         xi_p[1] = gp_xi[p][1];
780:         ConstructQ12D_Ni(xi_p,Ni_p);

782:         gp_x = 0.0;
783:         gp_y = 0.0;
784:         for (n = 0; n < NODES_PER_EL; n++) {
785:           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
786:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
787:         }
788:         element_props[j][i].gp_coords[2*p]   = gp_x;
789:         element_props[j][i].gp_coords[2*p+1] = gp_y;
790:       }
791:     }
792:   }

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

797:   for (j = sj; j < sj+ny; j++) {
798:     for (i = si; i < si+nx; i++) {
799:       PetscScalar              centroid_x = _prop_coords[j][i].x; /* centroids of cell */
800:       PetscScalar              centroid_y = _prop_coords[j][i].y;
801:       PETSC_UNUSED PetscScalar coord_x,coord_y;

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

806:         opts_E  = 1.0;
807:         opts_nu = 0.33;
808:         PetscOptionsGetScalar(NULL,NULL,"-iso_E",&opts_E,&flg);
809:         PetscOptionsGetScalar(NULL,NULL,"-iso_nu",&opts_nu,&flg);

811:         for (p = 0; p < GAUSS_POINTS; p++) {
812:           element_props[j][i].E[p]  = opts_E;
813:           element_props[j][i].nu[p] = opts_nu;

815:           element_props[j][i].fx[p] = 0.0;
816:           element_props[j][i].fy[p] = 0.0;
817:         }
818:       } else if (coefficient_structure == 1) { /* step */
819:         PetscScalar opts_E0,opts_nu0,opts_xc;
820:         PetscScalar opts_E1,opts_nu1;

822:         opts_E0  = opts_E1  = 1.0;
823:         opts_nu0 = opts_nu1 = 0.333;
824:         opts_xc  = 0.5;
825:         PetscOptionsGetScalar(NULL,NULL,"-step_E0",&opts_E0,&flg);
826:         PetscOptionsGetScalar(NULL,NULL,"-step_nu0",&opts_nu0,&flg);
827:         PetscOptionsGetScalar(NULL,NULL,"-step_E1",&opts_E1,&flg);
828:         PetscOptionsGetScalar(NULL,NULL,"-step_nu1",&opts_nu1,&flg);
829:         PetscOptionsGetScalar(NULL,NULL,"-step_xc",&opts_xc,&flg);

831:         for (p = 0; p < GAUSS_POINTS; p++) {
832:           coord_x = centroid_x;
833:           coord_y = centroid_y;
834:           if (use_gp_coords) {
835:             coord_x = element_props[j][i].gp_coords[2*p];
836:             coord_y = element_props[j][i].gp_coords[2*p+1];
837:           }

839:           element_props[j][i].E[p]  = opts_E0;
840:           element_props[j][i].nu[p] = opts_nu0;
841:           if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
842:             element_props[j][i].E[p]  = opts_E1;
843:             element_props[j][i].nu[p] = opts_nu1;
844:           }

846:           element_props[j][i].fx[p] = 0.0;
847:           element_props[j][i].fy[p] = 0.0;
848:         }
849:       } else if (coefficient_structure == 2) { /* brick */
850:         PetscReal values_E[10];
851:         PetscReal values_nu[10];
852:         PetscInt  nbricks,maxnbricks;
853:         PetscInt  index,span;
854:         PetscInt  jj;

856:         flg        = PETSC_FALSE;
857:         maxnbricks = 10;
858:         PetscOptionsGetRealArray(NULL,NULL, "-brick_E",values_E,&maxnbricks,&flg);
859:         nbricks    = maxnbricks;

862:         flg        = PETSC_FALSE;
863:         maxnbricks = 10;
864:         PetscOptionsGetRealArray(NULL,NULL, "-brick_nu",values_nu,&maxnbricks,&flg);

868:         span = 1;
869:         PetscOptionsGetInt(NULL,NULL,"-brick_span",&span,&flg);

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

876:         for (p = 0; p < GAUSS_POINTS; p++) {
877:           element_props[j][i].E[p]  = values_E[index];
878:           element_props[j][i].nu[p] = values_nu[index];
879:         }
880:       } else if (coefficient_structure == 3) { /* sponge */
881:         PetscScalar opts_E0,opts_nu0;
882:         PetscScalar opts_E1,opts_nu1;
883:         PetscInt    opts_t,opts_w;
884:         PetscInt    ii,jj,ci,cj;

886:         opts_E0  = opts_E1  = 1.0;
887:         opts_nu0 = opts_nu1 = 0.333;
888:         PetscOptionsGetScalar(NULL,NULL,"-sponge_E0",&opts_E0,&flg);
889:         PetscOptionsGetScalar(NULL,NULL,"-sponge_nu0",&opts_nu0,&flg);
890:         PetscOptionsGetScalar(NULL,NULL,"-sponge_E1",&opts_E1,&flg);
891:         PetscOptionsGetScalar(NULL,NULL,"-sponge_nu1",&opts_nu1,&flg);

893:         opts_t = opts_w = 1;
894:         PetscOptionsGetInt(NULL,NULL,"-sponge_t",&opts_t,&flg);
895:         PetscOptionsGetInt(NULL,NULL,"-sponge_w",&opts_w,&flg);

897:         ii = (i)/(opts_t+opts_w+opts_t);
898:         jj = (j)/(opts_t+opts_w+opts_t);

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

903:         for (p = 0; p < GAUSS_POINTS; p++) {
904:           element_props[j][i].E[p]  = opts_E0;
905:           element_props[j][i].nu[p] = opts_nu0;
906:         }
907:         if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
908:           if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
909:             for (p = 0; p < GAUSS_POINTS; p++) {
910:               element_props[j][i].E[p]  = opts_E1;
911:               element_props[j][i].nu[p] = opts_nu1;
912:             }
913:           }
914:         }
915:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
916:     }
917:   }
918:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

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

922:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
923:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
924:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);

926:   PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
927:   if (!no_view) {
928:     DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
929:     DMDACoordViewGnuplot2d(elas_da,"mesh");
930:   }

932:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
933:   DMCreateMatrix(elas_da,&A);
934:   DMGetCoordinates(elas_da,&vel_coords);
935:   MatNullSpaceCreateRigidBody(vel_coords,&matnull);
936:   MatSetNearNullSpace(A,matnull);
937:   MatNullSpaceDestroy(&matnull);
938:   MatCreateVecs(A,&f,&X);

940:   /* assemble A11 */
941:   MatZeroEntries(A);
942:   VecZeroEntries(f);

944:   AssembleA_Elasticity(A,elas_da,da_prop,properties);
945:   /* build force vector */
946:   AssembleF_Elasticity(f,elas_da,da_prop,properties);

948:   KSPCreate(PETSC_COMM_WORLD,&ksp_E);
949:   KSPSetOptionsPrefix(ksp_E,"elas_");  /* elasticity */

951:   PetscOptionsGetBool(NULL,NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
952:   /* solve */
953:   if (!use_nonsymbc) {
954:     Mat        AA;
955:     Vec        ff,XX;
956:     IS         is;
957:     VecScatter scat;

959:     DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
960:     VecDuplicate(ff,&XX);

962:     KSPSetOperators(ksp_E,AA,AA);
963:     KSPSetFromOptions(ksp_E);

965:     KSPSolve(ksp_E,ff,XX);

967:     /* push XX back into X */
968:     DMDABCApplyCompression(elas_da,NULL,X);

970:     VecScatterCreate(XX,NULL,X,is,&scat);
971:     VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
972:     VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
973:     VecScatterDestroy(&scat);

975:     MatDestroy(&AA);
976:     VecDestroy(&ff);
977:     VecDestroy(&XX);
978:     ISDestroy(&is);
979:   } else {
980:     DMDABCApplyCompression(elas_da,A,f);

982:     KSPSetOperators(ksp_E,A,A);
983:     KSPSetFromOptions(ksp_E);

985:     KSPSolve(ksp_E,f,X);
986:   }

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

991:   VecDestroy(&X);
992:   VecDestroy(&f);
993:   MatDestroy(&A);

995:   DMDestroy(&elas_da);
996:   DMDestroy(&da_prop);

998:   VecDestroy(&properties);
999:   VecDestroy(&l_properties);
1000:   return 0;
1001: }

1003: int main(int argc,char **args)
1004: {
1005:   PetscInt       mx,my;

1007:   PetscInitialize(&argc,&args,(char*)0,help);
1008:   mx   = my = 10;
1009:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1010:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1011:   solve_elasticity_2d(mx,my);
1012:   PetscFinalize();
1013:   return 0;
1014: }

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

1018: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1019: {
1020:   DM                     cda;
1021:   Vec                    coords;
1022:   PetscInt               si,sj,nx,ny,i,j;
1023:   PetscInt               M,N;
1024:   DMDACoor2d             **_coords;
1025:   const PetscInt         *g_idx;
1026:   PetscInt               *bc_global_ids;
1027:   PetscScalar            *bc_vals;
1028:   PetscInt               nbcs;
1029:   PetscInt               n_dofs;
1030:   ISLocalToGlobalMapping ltogm;

1033:   /* enforce bc's */
1034:   DMGetLocalToGlobalMapping(da,&ltogm);
1035:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1037:   DMGetCoordinateDM(da,&cda);
1038:   DMGetCoordinatesLocal(da,&coords);
1039:   DMDAVecGetArray(cda,coords,&_coords);
1040:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1041:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1043:   /* --- */

1045:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1046:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

1051:   i = nx-1;
1052:   for (j = 0; j < ny; j++) {
1053:     PetscInt                 local_id;
1054:     PETSC_UNUSED PetscScalar coordx,coordy;

1056:     local_id = i+j*nx;

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

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

1063:     bc_vals[j] =  bc_val;
1064:   }
1065:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1066:   nbcs = 0;
1067:   if ((si+nx) == (M)) nbcs = ny;

1069:   if (b) {
1070:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1071:     VecAssemblyBegin(b);
1072:     VecAssemblyEnd(b);
1073:   }
1074:   if (A) {
1075:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1076:   }

1078:   PetscFree(bc_vals);
1079:   PetscFree(bc_global_ids);

1081:   DMDAVecRestoreArray(cda,coords,&_coords);
1082:   return 0;
1083: }

1085: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1086: {
1087:   DM                     cda;
1088:   Vec                    coords;
1089:   PetscInt               si,sj,nx,ny,i,j;
1090:   PetscInt               M,N;
1091:   DMDACoor2d             **_coords;
1092:   const PetscInt         *g_idx;
1093:   PetscInt               *bc_global_ids;
1094:   PetscScalar            *bc_vals;
1095:   PetscInt               nbcs;
1096:   PetscInt               n_dofs;
1097:   ISLocalToGlobalMapping ltogm;

1100:   /* enforce bc's */
1101:   DMGetLocalToGlobalMapping(da,&ltogm);
1102:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1104:   DMGetCoordinateDM(da,&cda);
1105:   DMGetCoordinatesLocal(da,&coords);
1106:   DMDAVecGetArray(cda,coords,&_coords);
1107:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1108:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1110:   /* --- */

1112:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1113:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

1118:   i = 0;
1119:   for (j = 0; j < ny; j++) {
1120:     PetscInt                 local_id;
1121:     PETSC_UNUSED PetscScalar coordx,coordy;

1123:     local_id = i+j*nx;

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

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

1130:     bc_vals[j] =  bc_val;
1131:   }
1132:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1133:   nbcs = 0;
1134:   if (si == 0) nbcs = ny;

1136:   if (b) {
1137:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1138:     VecAssemblyBegin(b);
1139:     VecAssemblyEnd(b);
1140:   }
1141:   if (A) {
1142:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1143:   }

1145:   PetscFree(bc_vals);
1146:   PetscFree(bc_global_ids);

1148:   DMDAVecRestoreArray(cda,coords,&_coords);
1149:   return 0;
1150: }

1152: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1153: {
1155:   BCApply_EAST(elas_da,0,-1.0,A,f);
1156:   BCApply_EAST(elas_da,1, 0.0,A,f);
1157:   BCApply_WEST(elas_da,0,1.0,A,f);
1158:   BCApply_WEST(elas_da,1,0.0,A,f);
1159:   return 0;
1160: }

1162: static PetscErrorCode Orthogonalize(PetscInt n,Vec *vecs)
1163: {
1164:   PetscInt       i,j;
1165:   PetscScalar    dot;

1167:   for (i=0; i<n; i++) {
1168:   VecNormalize(vecs[i],NULL);
1169:      for (j=i+1; j<n; j++) {
1170:        VecDot(vecs[i],vecs[j],&dot);
1171:        VecAXPY(vecs[j],-dot,vecs[i]);
1172:      }
1173:   }
1174:   return 0;
1175: }

1177: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1178: {
1179:   PetscInt       start,end,m;
1180:   PetscInt       *unconstrained;
1181:   PetscInt       cnt,i;
1182:   Vec            x;
1183:   PetscScalar    *_x;
1184:   IS             is;
1185:   VecScatter     scat;

1188:   /* push bc's into f and A */
1189:   VecDuplicate(f,&x);
1190:   BCApply_EAST(elas_da,0,-1.0,A,x);
1191:   BCApply_EAST(elas_da,1, 0.0,A,x);
1192:   BCApply_WEST(elas_da,0,1.0,A,x);
1193:   BCApply_WEST(elas_da,1,0.0,A,x);

1195:   /* define which dofs are not constrained */
1196:   VecGetLocalSize(x,&m);
1197:   PetscMalloc1(m,&unconstrained);
1198:   VecGetOwnershipRange(x,&start,&end);
1199:   VecGetArray(x,&_x);
1200:   cnt  = 0;
1201:   for (i = 0; i < m; i+=2) {
1202:     PetscReal val1,val2;

1204:     val1 = PetscRealPart(_x[i]);
1205:     val2 = PetscRealPart(_x[i+1]);
1206:     if (PetscAbs(val1) < 0.1 && PetscAbs(val2) < 0.1) {
1207:       unconstrained[cnt] = start + i;
1208:       cnt++;
1209:       unconstrained[cnt] = start + i + 1;
1210:       cnt++;
1211:     }
1212:   }
1213:   VecRestoreArray(x,&_x);

1215:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
1216:   PetscFree(unconstrained);
1217:   ISSetBlockSize(is,2);

1219:   /* define correction for dirichlet in the rhs */
1220:   MatMult(A,x,f);
1221:   VecScale(f,-1.0);

1223:   /* get new matrix */
1224:   MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
1225:   /* get new vector */
1226:   MatCreateVecs(*AA,NULL,ff);

1228:   VecScatterCreate(f,is,*ff,NULL,&scat);
1229:   VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1230:   VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);

1232:   {                             /* Constrain near-null space */
1233:     PetscInt     nvecs;
1234:     const        Vec *vecs;
1235:     Vec          *uvecs;
1236:     PetscBool    has_const;
1237:     MatNullSpace mnull,unull;

1239:     MatGetNearNullSpace(A,&mnull);
1240:     MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);
1241:     VecDuplicateVecs(*ff,nvecs,&uvecs);
1242:     for (i=0; i<nvecs; i++) {
1243:       VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1244:       VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1245:     }
1246:     Orthogonalize(nvecs,uvecs);
1247:     MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,nvecs,uvecs,&unull);
1248:     MatSetNearNullSpace(*AA,unull);
1249:     MatNullSpaceDestroy(&unull);
1250:     VecDestroyVecs(nvecs,&uvecs);
1251:   }

1253:   VecScatterDestroy(&scat);

1255:   *dofs = is;
1256:   VecDestroy(&x);
1257:   return 0;
1258: }

1260: /*TEST

1262:    build:
1263:       requires: !complex !single

1265:    test:
1266:       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
1267:       output_file: output/ex49_1.out

1269:    test:
1270:       suffix: 2
1271:       nsize: 4
1272:       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

1274:    test:
1275:       suffix: 3
1276:       nsize: 4
1277:       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

1279:    test:
1280:       suffix: 4
1281:       nsize: 4
1282:       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
1283:       requires: ml

1285:    test:
1286:       suffix: 5
1287:       nsize: 3
1288:       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

1290:    test:
1291:       suffix: 6
1292:       nsize: 4
1293:       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

1295:    test:
1296:       suffix: 7
1297:       nsize: 4
1298:       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

1300:    test:
1301:       suffix: 8
1302:       nsize: 4
1303:       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

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

1310:    test:
1311:       nsize: 4
1312:       suffix: bddc
1313:       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

1315:    test:
1316:       nsize: 4
1317:       suffix: bddc_unsym
1318:       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

1320:    test:
1321:       nsize: 4
1322:       suffix: bddc_unsym_deluxe
1323:       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

1325:    test:
1326:       nsize: 4
1327:       suffix: fetidp_unsym_deluxe
1328:       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

1330:    test:
1331:       nsize: 4
1332:       suffix: bddc_layerjump
1333:       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

1335:    test:
1336:       nsize: 4
1337:       suffix: bddc_subdomainjump
1338:       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

1340:    test:
1341:       nsize: 9
1342:       suffix: bddc_subdomainjump_deluxe
1343:       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
1344: TEST*/