Actual source code: ex49.c

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

 64: /* Contributed by Dave May */

 66: #include <petscksp.h>
 67: #include <petscdm.h>
 68: #include <petscdmda.h>

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


 74: #define NSD            2 /* number of spatial dimensions */
 75: #define NODES_PER_EL   4 /* nodes per element */
 76: #define U_DOFS         2 /* degrees of freedom per displacement node */
 77: #define GAUSS_POINTS   4

 79: /* cell based evaluation */
 80: typedef struct {
 81:   PetscScalar E,nu,fx,fy;
 82: } Coefficients;

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

 93: typedef struct {
 94:   PetscScalar ux_dof;
 95:   PetscScalar uy_dof;
 96: } ElasticityDOF;


 99: /*

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

105:  B = [ d_dx   0   ]
106:      [  0    d_dy ]
107:      [ d_dy  d_dx ]

109:  */

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

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

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

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

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

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

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

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

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


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

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

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


192: /* procs to the left claim the ghost node as their element */
195: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
196: {
198:   PetscInt       m,n,p,M,N,P;
199:   PetscInt       sx,sy,sz;

202:   DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
203:   DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);

205:   if (mxl != NULL) {
206:     *mxl = m;
207:     if ((sx+m) == M) *mxl = m-1;    /* last proc */
208:   }
209:   if (myl != NULL) {
210:     *myl = n;
211:     if ((sy+n) == N) *myl = n-1;  /* last proc */
212:   }
213:   if (mzl != NULL) {
214:     *mzl = p;
215:     if ((sz+p) == P) *mzl = p-1;  /* last proc */
216:   }
217:   return(0);
218: }

222: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
223: {
225:   PetscInt       si,sj,sk;

228:   DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);

230:   if (sx) {
231:     *sx = si;
232:     if (si != 0) *sx = si+1;
233:   }
234:   if (sy) {
235:     *sy = sj;
236:     if (sj != 0) *sy = sj+1;
237:   }

239:   if (sk) {
240:     *sz = sk;
241:     if (sk != 0) *sz = sk+1;
242:   }

244:   DMDAGetLocalElementSize(da,mx,my,mz);
245:   return(0);
246: }

250: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
251: {
253:   PetscMPIInt    rank;
254:   PetscInt       proc_I,proc_J;
255:   PetscInt       cpu_x,cpu_y;
256:   PetscInt       local_mx,local_my;
257:   Vec            vlx,vly;
258:   PetscInt       *LX,*LY,i;
259:   PetscScalar    *_a;
260:   Vec            V_SEQ;
261:   VecScatter     ctx;

264:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

268:   proc_J = rank/cpu_x;
269:   proc_I = rank-cpu_x*proc_J;

271:   PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
272:   PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);

274:   DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
275:   VecCreate(PETSC_COMM_WORLD,&vlx);
276:   VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
277:   VecSetFromOptions(vlx);

279:   VecCreate(PETSC_COMM_WORLD,&vly);
280:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
281:   VecSetFromOptions(vly);

283:   VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
284:   VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
285:   VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
286:   VecAssemblyBegin(vly);VecAssemblyEnd(vly);



290:   VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
291:   VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
292:   VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
293:   VecGetArray(V_SEQ,&_a);
294:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
295:   VecRestoreArray(V_SEQ,&_a);
296:   VecScatterDestroy(&ctx);
297:   VecDestroy(&V_SEQ);

299:   VecScatterCreateToAll(vly,&ctx,&V_SEQ);
300:   VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
301:   VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
302:   VecGetArray(V_SEQ,&_a);
303:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
304:   VecRestoreArray(V_SEQ,&_a);
305:   VecScatterDestroy(&ctx);
306:   VecDestroy(&V_SEQ);

308:   *_lx = LX;
309:   *_ly = LY;

311:   VecDestroy(&vlx);
312:   VecDestroy(&vly);
313:   return(0);
314: }

318: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
319: {
320:   DM             cda;
321:   Vec            coords;
322:   DMDACoor2d     **_coords;
323:   PetscInt       si,sj,nx,ny,i,j;
324:   FILE           *fp;
325:   char           fname[PETSC_MAX_PATH_LEN];
326:   PetscMPIInt    rank;

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

336:   DMGetCoordinateDM(da,&cda);
337:   DMGetCoordinatesLocal(da,&coords);
338:   DMDAVecGetArray(cda,coords,&_coords);
339:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
340:   for (j = sj; j < sj+ny-1; j++) {
341:     for (i = si; i < si+nx-1; i++) {
342:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
343:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
344:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
345:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
346:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
347:     }
348:   }
349:   DMDAVecRestoreArray(cda,coords,&_coords);

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

357: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
358: {
359:   DM             cda;
360:   Vec            coords,local_fields;
361:   DMDACoor2d     **_coords;
362:   FILE           *fp;
363:   char           fname[PETSC_MAX_PATH_LEN];
364:   const char     *field_name;
365:   PetscMPIInt    rank;
366:   PetscInt       si,sj,nx,ny,i,j;
367:   PetscInt       n_dofs,d;
368:   PetscScalar    *_fields;

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

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


387:   DMGetCoordinateDM(da,&cda);
388:   DMGetCoordinatesLocal(da,&coords);
389:   DMDAVecGetArray(cda,coords,&_coords);
390:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

392:   DMCreateLocalVector(da,&local_fields);
393:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
394:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
395:   VecGetArray(local_fields,&_fields);

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

402:       coord_x = _coords[j][i].x;
403:       coord_y = _coords[j][i].y;

405:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
406:       for (d = 0; d < n_dofs; d++) {
407:         field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
408:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
409:       }
410:       PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
411:     }
412:   }
413:   VecRestoreArray(local_fields,&_fields);
414:   VecDestroy(&local_fields);

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

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

424: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
425: {
426:   DM                     cda;
427:   Vec                    local_fields;
428:   FILE                   *fp;
429:   char                   fname[PETSC_MAX_PATH_LEN];
430:   const char             *field_name;
431:   PetscMPIInt            rank;
432:   PetscInt               si,sj,nx,ny,i,j,p;
433:   PetscInt               n_dofs,d;
434:   GaussPointCoefficients **_coefficients;
435:   PetscErrorCode         ierr;

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

443:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
444:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
445:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
446:   for (d = 0; d < n_dofs; d++) {
447:     DMDAGetFieldName(da,d,&field_name);
448:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
449:   }
450:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");


453:   DMGetCoordinateDM(da,&cda);
454:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

456:   DMCreateLocalVector(da,&local_fields);
457:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
458:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
459:   DMDAVecGetArray(da,local_fields,&_coefficients);

461:   for (j = sj; j < sj+ny; j++) {
462:     for (i = si; i < si+nx; i++) {
463:       PetscScalar coord_x,coord_y;

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

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

471:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e %1.6e",
472:                             PetscRealPart(_coefficients[j][i].E[p]),PetscRealPart(_coefficients[j][i].nu[p]),
473:                             PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
474:         PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
475:       }
476:     }
477:   }
478:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
479:   VecDestroy(&local_fields);

481:   PetscFClose(PETSC_COMM_SELF,fp);
482:   return(0);
483: }

485: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar E[],PetscScalar nu[])
486: {
487:   PetscInt    ngp;
488:   PetscScalar gp_xi[GAUSS_POINTS][2];
489:   PetscScalar gp_weight[GAUSS_POINTS];
490:   PetscInt    p,i,j,k,l;
491:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
492:   PetscScalar J_p;
493:   PetscScalar B[3][U_DOFS*NODES_PER_EL];
494:   PetscScalar prop_E,prop_nu,factor,constit_D[3][3];

496:   /* define quadrature rule */
497:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

499:   /* evaluate integral */
500:   for (p = 0; p < ngp; p++) {
501:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
502:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

504:     for (i = 0; i < NODES_PER_EL; i++) {
505:       PetscScalar d_dx_i = GNx_p[0][i];
506:       PetscScalar d_dy_i = GNx_p[1][i];

508:       B[0][2*i] = d_dx_i;  B[0][2*i+1] = 0.0;
509:       B[1][2*i] = 0.0;     B[1][2*i+1] = d_dy_i;
510:       B[2][2*i] = d_dy_i;  B[2][2*i+1] = d_dx_i;
511:     }

513:     /* form D for the quadrature point */
514:     prop_E          = E[p];
515:     prop_nu         = nu[p];
516:     factor          = prop_E / ((1.0+prop_nu)*(1.0-2.0*prop_nu));
517:     constit_D[0][0] = 1.0-prop_nu;  constit_D[0][1] = prop_nu;      constit_D[0][2] = 0.0;
518:     constit_D[1][0] = prop_nu;      constit_D[1][1] = 1.0-prop_nu;  constit_D[1][2] = 0.0;
519:     constit_D[2][0] = 0.0;          constit_D[2][1] = 0.0;          constit_D[2][2] = 0.5*(1.0-2.0*prop_nu);
520:     for (i = 0; i < 3; i++) {
521:       for (j = 0; j < 3; j++) {
522:         constit_D[i][j] = factor * constit_D[i][j] * gp_weight[p] * J_p;
523:       }
524:     }

526:     /* form Bt tildeD B */
527:     /*
528:      Ke_ij = Bt_ik . D_kl . B_lj
529:      = B_ki . D_kl . B_lj
530:      */
531:     for (i = 0; i < 8; i++) {
532:       for (j = 0; j < 8; j++) {
533:         for (k = 0; k < 3; k++) {
534:           for (l = 0; l < 3; l++) {
535:             Ke[8*i+j] = Ke[8*i+j] + B[k][i] * constit_D[k][l] * B[l][j];
536:           }
537:         }
538:       }
539:     }

541:   } /* end quadrature */
542: }

544: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
545: {
546:   PetscInt    ngp;
547:   PetscScalar gp_xi[GAUSS_POINTS][2];
548:   PetscScalar gp_weight[GAUSS_POINTS];
549:   PetscInt    p,i;
550:   PetscScalar Ni_p[NODES_PER_EL];
551:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
552:   PetscScalar J_p,fac;

554:   /* define quadrature rule */
555:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

557:   /* evaluate integral */
558:   for (p = 0; p < ngp; p++) {
559:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
560:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
561:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
562:     fac = gp_weight[p]*J_p;

564:     for (i = 0; i < NODES_PER_EL; i++) {
565:       Fe[NSD*i]   += fac*Ni_p[i]*fx[p];
566:       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
567:     }
568:   }
569: }

571: /*
572:  i,j are the element indices
573:  The unknown is a vector quantity.
574:  The s[].c is used to indicate the degree of freedom.
575:  */
578: static PetscErrorCode DMDAGetElementEqnums_u(MatStencil s_u[],PetscInt i,PetscInt j)
579: {
581:   /* displacement */
582:   /* node 0 */
583:   s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0;          /* Ux0 */
584:   s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1;          /* Uy0 */

586:   /* node 1 */
587:   s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0;        /* Ux1 */
588:   s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1;        /* Uy1 */

590:   /* node 2 */
591:   s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0;      /* Ux2 */
592:   s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1;      /* Uy2 */

594:   /* node 3 */
595:   s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0;        /* Ux3 */
596:   s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1;        /* Uy3 */
597:   return(0);
598: }

602: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
603: {
605:   /* get coords for the element */
606:   el_coords[NSD*0+0] = _coords[ej][ei].x;      el_coords[NSD*0+1] = _coords[ej][ei].y;
607:   el_coords[NSD*1+0] = _coords[ej+1][ei].x;    el_coords[NSD*1+1] = _coords[ej+1][ei].y;
608:   el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;  el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
609:   el_coords[NSD*3+0] = _coords[ej][ei+1].x;    el_coords[NSD*3+1] = _coords[ej][ei+1].y;
610:   return(0);
611: }

615: static PetscErrorCode AssembleA_Elasticity(Mat A,DM elas_da,DM properties_da,Vec properties)
616: {
617:   DM                     cda;
618:   Vec                    coords;
619:   DMDACoor2d             **_coords;
620:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
621:   PetscInt               sex,sey,mx,my;
622:   PetscInt               ei,ej;
623:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
624:   PetscScalar            el_coords[NODES_PER_EL*NSD];
625:   Vec                    local_properties;
626:   GaussPointCoefficients **props;
627:   PetscScalar            *prop_E,*prop_nu;
628:   PetscErrorCode         ierr;

631:   /* setup for coords */
632:   DMGetCoordinateDM(elas_da,&cda);
633:   DMGetCoordinatesLocal(elas_da,&coords);
634:   DMDAVecGetArray(cda,coords,&_coords);

636:   /* setup for coefficients */
637:   DMCreateLocalVector(properties_da,&local_properties);
638:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
639:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
640:   DMDAVecGetArray(properties_da,local_properties,&props);

642:   DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
643:   for (ej = sey; ej < sey+my; ej++) {
644:     for (ei = sex; ei < sex+mx; ei++) {
645:       /* get coords for the element */
646:       GetElementCoords(_coords,ei,ej,el_coords);

648:       /* get coefficients for the element */
649:       prop_E  = props[ej][ei].E;
650:       prop_nu = props[ej][ei].nu;

652:       /* initialise element stiffness matrix */
653:       PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);

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

658:       /* insert element matrix into global matrix */
659:       DMDAGetElementEqnums_u(u_eqn,ei,ej);
660:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
661:     }
662:   }
663:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
664:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

668:   DMDAVecRestoreArray(properties_da,local_properties,&props);
669:   VecDestroy(&local_properties);
670:   return(0);
671: }


676: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F,MatStencil u_eqn[],PetscScalar Fe_u[])
677: {
678:   PetscInt n;

681:   for (n = 0; n < 4; n++) {
682:     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];
683:     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];
684:   }
685:   return(0);
686: }

690: static PetscErrorCode AssembleF_Elasticity(Vec F,DM elas_da,DM properties_da,Vec properties)
691: {
692:   DM                     cda;
693:   Vec                    coords;
694:   DMDACoor2d             **_coords;
695:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
696:   PetscInt               sex,sey,mx,my;
697:   PetscInt               ei,ej;
698:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
699:   PetscScalar            el_coords[NODES_PER_EL*NSD];
700:   Vec                    local_properties;
701:   GaussPointCoefficients **props;
702:   PetscScalar            *prop_fx,*prop_fy;
703:   Vec                    local_F;
704:   ElasticityDOF          **ff;
705:   PetscErrorCode         ierr;

708:   /* setup for coords */
709:   DMGetCoordinateDM(elas_da,&cda);
710:   DMGetCoordinatesLocal(elas_da,&coords);
711:   DMDAVecGetArray(cda,coords,&_coords);

713:   /* setup for coefficients */
714:   DMGetLocalVector(properties_da,&local_properties);
715:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
716:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
717:   DMDAVecGetArray(properties_da,local_properties,&props);

719:   /* get acces to the vector */
720:   DMGetLocalVector(elas_da,&local_F);
721:   VecZeroEntries(local_F);
722:   DMDAVecGetArray(elas_da,local_F,&ff);


725:   DMDAGetElementCorners(elas_da,&sex,&sey,0,&mx,&my,0);
726:   for (ej = sey; ej < sey+my; ej++) {
727:     for (ei = sex; ei < sex+mx; ei++) {
728:       /* get coords for the element */
729:       GetElementCoords(_coords,ei,ej,el_coords);

731:       /* get coefficients for the element */
732:       prop_fx = props[ej][ei].fx;
733:       prop_fy = props[ej][ei].fy;

735:       /* initialise element stiffness matrix */
736:       PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);

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

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

744:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,Fe);
745:     }
746:   }

748:   DMDAVecRestoreArray(elas_da,local_F,&ff);
749:   DMLocalToGlobalBegin(elas_da,local_F,ADD_VALUES,F);
750:   DMLocalToGlobalEnd(elas_da,local_F,ADD_VALUES,F);
751:   DMRestoreLocalVector(elas_da,&local_F);

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

755:   DMDAVecRestoreArray(properties_da,local_properties,&props);
756:   DMRestoreLocalVector(properties_da,&local_properties);
757:   return(0);
758: }

762: static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
763: {
764:   DM                     elas_da,da_prop;
765:   PetscInt               u_dof,dof,stencil_width;
766:   Mat                    A;
767:   PetscInt               mxl,myl;
768:   DM                     prop_cda,vel_cda;
769:   Vec                    prop_coords,vel_coords;
770:   PetscInt               si,sj,nx,ny,i,j,p;
771:   Vec                    f,X;
772:   PetscInt               prop_dof,prop_stencil_width;
773:   Vec                    properties,l_properties;
774:   MatNullSpace           matnull;
775:   PetscReal              dx,dy;
776:   PetscInt               M,N;
777:   DMDACoor2d             **_prop_coords,**_vel_coords;
778:   GaussPointCoefficients **element_props;
779:   KSP                    ksp_E;
780:   PetscInt               coefficient_structure = 0;
781:   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
782:   PetscBool              use_gp_coords = PETSC_FALSE;
783:   PetscBool              use_nonsymbc  = PETSC_FALSE;
784:   PetscBool              no_view       = PETSC_FALSE;
785:   PetscBool              flg;
786:   PetscErrorCode         ierr;

789:   /* Generate the da for velocity and pressure */
790:   /*
791:    We use Q1 elements for the temperature.
792:    FEM has a 9-point stencil (BOX) or connectivity pattern
793:    Num nodes in each direction is mx+1, my+1
794:    */
795:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
796:   dof           = u_dof;
797:   stencil_width = 1;
798:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
799:                                mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&elas_da);
800:   DMDASetFieldName(elas_da,0,"Ux");
801:   DMDASetFieldName(elas_da,1,"Uy");

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


807:   /* Generate element properties, we will assume all material properties are constant over the element */
808:   /* local number of elements */
809:   DMDAGetLocalElementSize(elas_da,&mxl,&myl,NULL);

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

815:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
816:   prop_stencil_width = 0;
817:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
818:                                     mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
819:   PetscFree(lx);
820:   PetscFree(ly);

822:   /* define centroid positions */
823:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
824:   dx   = 1.0/((PetscReal)(M));
825:   dy   = 1.0/((PetscReal)(N));

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

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

832:   DMCreateGlobalVector(da_prop,&properties);
833:   DMCreateLocalVector(da_prop,&l_properties);
834:   DMDAVecGetArray(da_prop,l_properties,&element_props);

836:   DMGetCoordinateDM(da_prop,&prop_cda);
837:   DMGetCoordinatesLocal(da_prop,&prop_coords);
838:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

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

842:   DMGetCoordinateDM(elas_da,&vel_cda);
843:   DMGetCoordinatesLocal(elas_da,&vel_coords);
844:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);


847:   /* interpolate the coordinates */
848:   for (j = sj; j < sj+ny; j++) {
849:     for (i = si; i < si+nx; i++) {
850:       PetscInt    ngp;
851:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
852:       PetscScalar el_coords[8];

854:       GetElementCoords(_vel_coords,i,j,el_coords);
855:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

857:       for (p = 0; p < GAUSS_POINTS; p++) {
858:         PetscScalar gp_x,gp_y;
859:         PetscInt    n;
860:         PetscScalar xi_p[2],Ni_p[4];

862:         xi_p[0] = gp_xi[p][0];
863:         xi_p[1] = gp_xi[p][1];
864:         ConstructQ12D_Ni(xi_p,Ni_p);

866:         gp_x = 0.0;
867:         gp_y = 0.0;
868:         for (n = 0; n < NODES_PER_EL; n++) {
869:           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
870:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
871:         }
872:         element_props[j][i].gp_coords[2*p]   = gp_x;
873:         element_props[j][i].gp_coords[2*p+1] = gp_y;
874:       }
875:     }
876:   }

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

881:   for (j = sj; j < sj+ny; j++) {
882:     for (i = si; i < si+nx; i++) {
883:       PetscScalar              centroid_x = _prop_coords[j][i].x; /* centroids of cell */
884:       PetscScalar              centroid_y = _prop_coords[j][i].y;
885:       PETSC_UNUSED PetscScalar coord_x,coord_y;


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

891:         opts_E  = 1.0;
892:         opts_nu = 0.33;
893:         PetscOptionsGetScalar(NULL,"-iso_E",&opts_E,&flg);
894:         PetscOptionsGetScalar(NULL,"-iso_nu",&opts_nu,&flg);

896:         for (p = 0; p < GAUSS_POINTS; p++) {
897:           element_props[j][i].E[p]  = opts_E;
898:           element_props[j][i].nu[p] = opts_nu;

900:           element_props[j][i].fx[p] = 0.0;
901:           element_props[j][i].fy[p] = 0.0;
902:         }
903:       } else if (coefficient_structure == 1) { /* step */
904:         PetscScalar opts_E0,opts_nu0,opts_xc;
905:         PetscScalar opts_E1,opts_nu1;

907:         opts_E0  = opts_E1  = 1.0;
908:         opts_nu0 = opts_nu1 = 0.333;
909:         opts_xc  = 0.5;
910:         PetscOptionsGetScalar(NULL,"-step_E0",&opts_E0,&flg);
911:         PetscOptionsGetScalar(NULL,"-step_nu0",&opts_nu0,&flg);
912:         PetscOptionsGetScalar(NULL,"-step_E1",&opts_E1,&flg);
913:         PetscOptionsGetScalar(NULL,"-step_nu1",&opts_nu1,&flg);
914:         PetscOptionsGetScalar(NULL,"-step_xc",&opts_xc,&flg);

916:         for (p = 0; p < GAUSS_POINTS; p++) {
917:           coord_x = centroid_x;
918:           coord_y = centroid_y;
919:           if (use_gp_coords) {
920:             coord_x = element_props[j][i].gp_coords[2*p];
921:             coord_y = element_props[j][i].gp_coords[2*p+1];
922:           }

924:           element_props[j][i].E[p]  = opts_E0;
925:           element_props[j][i].nu[p] = opts_nu0;
926:           if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
927:             element_props[j][i].E[p]  = opts_E1;
928:             element_props[j][i].nu[p] = opts_nu1;
929:           }

931:           element_props[j][i].fx[p] = 0.0;
932:           element_props[j][i].fy[p] = 0.0;
933:         }
934:       } else if (coefficient_structure == 2) { /* brick */
935:         PetscReal values_E[10];
936:         PetscReal values_nu[10];
937:         PetscInt  nbricks,maxnbricks;
938:         PetscInt  index,span;
939:         PetscInt  jj;

941:         flg        = PETSC_FALSE;
942:         maxnbricks = 10;
943:         PetscOptionsGetRealArray(NULL, "-brick_E",values_E,&maxnbricks,&flg);
944:         nbricks    = maxnbricks;
945:         if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of E values for each brick");

947:         flg        = PETSC_FALSE;
948:         maxnbricks = 10;
949:         PetscOptionsGetRealArray(NULL, "-brick_nu",values_nu,&maxnbricks,&flg);
950:         if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply a list of nu values for each brick");
951:         if (maxnbricks != nbricks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"User must supply equal numbers of values for E and nu");

953:         span = 1;
954:         PetscOptionsGetInt(NULL,"-brick_span",&span,&flg);

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

961:         for (p = 0; p < GAUSS_POINTS; p++) {
962:           element_props[j][i].E[p]  = values_E[index];
963:           element_props[j][i].nu[p] = values_nu[index];
964:         }
965:       } else if (coefficient_structure == 3) { /* sponge */
966:         PetscScalar opts_E0,opts_nu0;
967:         PetscScalar opts_E1,opts_nu1;
968:         PetscInt    opts_t,opts_w;
969:         PetscInt    ii,jj,ci,cj;

971:         opts_E0  = opts_E1  = 1.0;
972:         opts_nu0 = opts_nu1 = 0.333;
973:         PetscOptionsGetScalar(NULL,"-sponge_E0",&opts_E0,&flg);
974:         PetscOptionsGetScalar(NULL,"-sponge_nu0",&opts_nu0,&flg);
975:         PetscOptionsGetScalar(NULL,"-sponge_E1",&opts_E1,&flg);
976:         PetscOptionsGetScalar(NULL,"-sponge_nu1",&opts_nu1,&flg);

978:         opts_t = opts_w = 1;
979:         PetscOptionsGetInt(NULL,"-sponge_t",&opts_t,&flg);
980:         PetscOptionsGetInt(NULL,"-sponge_w",&opts_w,&flg);

982:         ii = (i)/(opts_t+opts_w+opts_t);
983:         jj = (j)/(opts_t+opts_w+opts_t);

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

988:         for (p = 0; p < GAUSS_POINTS; p++) {
989:           element_props[j][i].E[p]  = opts_E0;
990:           element_props[j][i].nu[p] = opts_nu0;
991:         }
992:         if ((ci >= opts_t) && (ci < opts_t+opts_w)) {
993:           if ((cj >= opts_t) && (cj < opts_t+opts_w)) {
994:             for (p = 0; p < GAUSS_POINTS; p++) {
995:               element_props[j][i].E[p]  = opts_E1;
996:               element_props[j][i].nu[p] = opts_nu1;
997:             }
998:           }
999:         }

1001:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1002:     }
1003:   }
1004:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

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

1008:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1009:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1010:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);

1012:   PetscOptionsGetBool(NULL,"-no_view",&no_view,NULL);
1013:   if (!no_view) {
1014:     DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");
1015:     DMDACoordViewGnuplot2d(elas_da,"mesh");
1016:   }

1018:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1019:   DMSetMatType(elas_da,MATAIJ);
1020:   DMCreateMatrix(elas_da,&A);
1021:   DMGetCoordinates(elas_da,&vel_coords);
1022:   MatNullSpaceCreateRigidBody(vel_coords,&matnull);
1023:   MatSetNearNullSpace(A,matnull);
1024:   MatNullSpaceDestroy(&matnull);
1025:   MatGetVecs(A,&f,&X);

1027:   /* assemble A11 */
1028:   MatZeroEntries(A);
1029:   VecZeroEntries(f);

1031:   AssembleA_Elasticity(A,elas_da,da_prop,properties);
1032:   /* build force vector */
1033:   AssembleF_Elasticity(f,elas_da,da_prop,properties);


1036:   KSPCreate(PETSC_COMM_WORLD,&ksp_E);
1037:   KSPSetOptionsPrefix(ksp_E,"elas_");  /* elasticity */

1039:   PetscOptionsGetBool(NULL,"-use_nonsymbc",&use_nonsymbc,&flg);
1040:   /* solve */
1041:   if (!use_nonsymbc) {
1042:     Mat        AA;
1043:     Vec        ff,XX;
1044:     IS         is;
1045:     VecScatter scat;

1047:     DMDABCApplySymmetricCompression(elas_da,A,f,&is,&AA,&ff);
1048:     VecDuplicate(ff,&XX);

1050:     KSPSetOperators(ksp_E,AA,AA);
1051:     KSPSetFromOptions(ksp_E);

1053:     KSPSolve(ksp_E,ff,XX);

1055:     /* push XX back into X */
1056:     DMDABCApplyCompression(elas_da,NULL,X);

1058:     VecScatterCreate(XX,NULL,X,is,&scat);
1059:     VecScatterBegin(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1060:     VecScatterEnd(scat,XX,X,INSERT_VALUES,SCATTER_FORWARD);
1061:     VecScatterDestroy(&scat);

1063:     MatDestroy(&AA);
1064:     VecDestroy(&ff);
1065:     VecDestroy(&XX);
1066:     ISDestroy(&is);
1067:   } else {
1068:     DMDABCApplyCompression(elas_da,A,f);

1070:     KSPSetOperators(ksp_E,A,A);
1071:     KSPSetFromOptions(ksp_E);

1073:     KSPSolve(ksp_E,f,X);
1074:   }

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

1079:   VecDestroy(&X);
1080:   VecDestroy(&f);
1081:   MatDestroy(&A);

1083:   DMDestroy(&elas_da);
1084:   DMDestroy(&da_prop);

1086:   VecDestroy(&properties);
1087:   VecDestroy(&l_properties);
1088:   return(0);
1089: }

1093: int main(int argc,char **args)
1094: {
1096:   PetscInt       mx,my;

1098:   PetscInitialize(&argc,&args,(char*)0,help);

1100:   mx   = my = 10;
1101:   PetscOptionsGetInt(NULL,"-mx",&mx,NULL);
1102:   PetscOptionsGetInt(NULL,"-my",&my,NULL);

1104:   solve_elasticity_2d(mx,my);

1106:   PetscFinalize();
1107:   return 0;
1108: }

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

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

1130:   /* enforce bc's */
1131:   DMGetLocalToGlobalMapping(da,&ltogm);
1132:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

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

1140:   /* --- */

1142:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1143:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

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

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

1153:     local_id = i+j*nx;

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

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

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

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

1175:   PetscFree(bc_vals);
1176:   PetscFree(bc_global_ids);

1178:   DMDAVecRestoreArray(cda,coords,&_coords);
1179:   return(0);
1180: }

1184: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1185: {
1186:   DM                     cda;
1187:   Vec                    coords;
1188:   PetscInt               si,sj,nx,ny,i,j;
1189:   PetscInt               M,N;
1190:   DMDACoor2d             **_coords;
1191:   const PetscInt         *g_idx;
1192:   PetscInt               *bc_global_ids;
1193:   PetscScalar            *bc_vals;
1194:   PetscInt               nbcs;
1195:   PetscInt               n_dofs;
1196:   PetscErrorCode         ierr;
1197:   ISLocalToGlobalMapping ltogm;

1200:   /* enforce bc's */
1201:   DMGetLocalToGlobalMapping(da,&ltogm);
1202:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1204:   DMGetCoordinateDM(da,&cda);
1205:   DMGetCoordinatesLocal(da,&coords);
1206:   DMDAVecGetArray(cda,coords,&_coords);
1207:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1208:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1210:   /* --- */

1212:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1213:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

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

1218:   i = 0;
1219:   for (j = 0; j < ny; j++) {
1220:     PetscInt                 local_id;
1221:     PETSC_UNUSED PetscScalar coordx,coordy;

1223:     local_id = i+j*nx;

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

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

1230:     bc_vals[j] =  bc_val;
1231:   }
1232:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1233:   nbcs = 0;
1234:   if (si == 0) nbcs = ny;

1236:   if (b) {
1237:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1238:     VecAssemblyBegin(b);
1239:     VecAssemblyEnd(b);
1240:   }
1241:   if (A) {
1242:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1243:   }

1245:   PetscFree(bc_vals);
1246:   PetscFree(bc_global_ids);

1248:   DMDAVecRestoreArray(cda,coords,&_coords);
1249:   return(0);
1250: }

1254: static PetscErrorCode DMDABCApplyCompression(DM elas_da,Mat A,Vec f)
1255: {

1259:   BCApply_EAST(elas_da,0,-1.0,A,f);
1260:   BCApply_EAST(elas_da,1, 0.0,A,f);
1261:   BCApply_WEST(elas_da,0,1.0,A,f);
1262:   BCApply_WEST(elas_da,1,0.0,A,f);
1263:   return(0);
1264: }

1268: static PetscErrorCode DMDABCApplySymmetricCompression(DM elas_da,Mat A,Vec f,IS *dofs,Mat *AA,Vec *ff)
1269: {
1271:   PetscInt       start,end,m;
1272:   PetscInt       *unconstrained;
1273:   PetscInt       cnt,i;
1274:   Vec            x;
1275:   PetscScalar    *_x;
1276:   IS             is;
1277:   VecScatter     scat;

1280:   /* push bc's into f and A */
1281:   VecDuplicate(f,&x);
1282:   BCApply_EAST(elas_da,0,-1.0,A,x);
1283:   BCApply_EAST(elas_da,1, 0.0,A,x);
1284:   BCApply_WEST(elas_da,0,1.0,A,x);
1285:   BCApply_WEST(elas_da,1,0.0,A,x);

1287:   /* define which dofs are not constrained */
1288:   VecGetLocalSize(x,&m);
1289:   PetscMalloc(sizeof(PetscInt)*m,&unconstrained);
1290:   VecGetOwnershipRange(x,&start,&end);
1291:   VecGetArray(x,&_x);
1292:   cnt  = 0;
1293:   for (i = 0; i < m; i++) {
1294:     PetscReal val;

1296:     val = PetscRealPart(_x[i]);
1297:     if (fabs(val) < 0.1) {
1298:       unconstrained[cnt] = start + i;
1299:       cnt++;
1300:     }
1301:   }
1302:   VecRestoreArray(x,&_x);

1304:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,unconstrained,PETSC_COPY_VALUES,&is);
1305:   PetscFree(unconstrained);

1307:   /* define correction for dirichlet in the rhs */
1308:   MatMult(A,x,f);
1309:   VecScale(f,-1.0);

1311:   /* get new matrix */
1312:   MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,AA);
1313:   /* get new vector */
1314:   MatGetVecs(*AA,NULL,ff);

1316:   VecScatterCreate(f,is,*ff,NULL,&scat);
1317:   VecScatterBegin(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);
1318:   VecScatterEnd(scat,f,*ff,INSERT_VALUES,SCATTER_FORWARD);

1320:   {                             /* Constrain near-null space */
1321:     PetscInt nvecs;
1322:     const Vec *vecs;
1323:     Vec *uvecs;
1324:     PetscBool has_const;
1325:     MatNullSpace mnull,unull;
1326:     MatGetNearNullSpace(A,&mnull);
1327:     MatNullSpaceGetVecs(mnull,&has_const,&nvecs,&vecs);
1328:     VecDuplicateVecs(*ff,nvecs,&uvecs);
1329:     for (i=0; i<nvecs; i++) {
1330:       VecScatterBegin(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1331:       VecScatterEnd(scat,vecs[i],uvecs[i],INSERT_VALUES,SCATTER_FORWARD);
1332:     }
1333:     MatNullSpaceCreate(PetscObjectComm((PetscObject)A),has_const,nvecs,uvecs,&unull);
1334:     MatSetNearNullSpace(*AA,unull);
1335:     MatNullSpaceDestroy(&unull);
1336:     for (i=0; i<nvecs; i++) {
1337:       VecDestroy(&uvecs[i]);
1338:     }
1339:     PetscFree(uvecs);
1340:   }

1342:   VecScatterDestroy(&scat);

1344:   *dofs = is;
1345:   VecDestroy(&x);
1346:   return(0);
1347: }