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 evaluated 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;
183:   gp_xi[0][1]  = -0.57735026919;
184:   gp_xi[1][0]  = -0.57735026919;
185:   gp_xi[1][1]  = 0.57735026919;
186:   gp_xi[2][0]  = 0.57735026919;
187:   gp_xi[2][1]  = 0.57735026919;
188:   gp_xi[3][0]  = 0.57735026919;
189:   gp_xi[3][1]  = -0.57735026919;
190:   gp_weight[0] = 1.0;
191:   gp_weight[1] = 1.0;
192:   gp_weight[2] = 1.0;
193:   gp_weight[3] = 1.0;
194: }

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

209:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

323:   DMGetCoordinateDM(da, &cda);
324:   DMGetCoordinatesLocal(da, &coords);
325:   DMDAVecGetArray(cda, coords, &_coords);
326:   DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);

328:   DMCreateLocalVector(da, &local_fields);
329:   DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
330:   DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
331:   VecGetArray(local_fields, &_fields);

333:   for (j = sj; j < sj + ny; j++) {
334:     for (i = si; i < si + nx; i++) {
335:       PetscScalar coord_x, coord_y;
336:       PetscScalar field_d;

338:       coord_x = _coords[j][i].x;
339:       coord_y = _coords[j][i].y;

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

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

354:   PetscFClose(PETSC_COMM_SELF, fp);
355:   return 0;
356: }

358: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
359: {
360:   DM                       cda;
361:   Vec                      local_fields;
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, p;
367:   PetscInt                 n_dofs, d;
368:   GaussPointCoefficients **_coefficients;

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

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

385:   DMGetCoordinateDM(da, &cda);
386:   DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);

388:   DMCreateLocalVector(da, &local_fields);
389:   DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
390:   DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
391:   DMDAVecGetArray(da, local_fields, &_coefficients);

393:   for (j = sj; j < sj + ny; j++) {
394:     for (i = si; i < si + nx; i++) {
395:       PetscScalar coord_x, coord_y;

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

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

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

411:   PetscFClose(PETSC_COMM_SELF, fp);
412:   return 0;
413: }

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

426:   /* define quadrature rule */
427:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

434:     for (i = 0; i < NODES_PER_EL; i++) {
435:       PetscScalar d_dx_i = GNx_p[0][i];
436:       PetscScalar d_dy_i = GNx_p[1][i];

438:       B[0][2 * i]     = d_dx_i;
439:       B[0][2 * i + 1] = 0.0;
440:       B[1][2 * i]     = 0.0;
441:       B[1][2 * i + 1] = d_dy_i;
442:       B[2][2 * i]     = d_dy_i;
443:       B[2][2 * i + 1] = d_dx_i;
444:     }

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

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

476:   } /* end quadrature */
477: }

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

489:   /* define quadrature rule */
490:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

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

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

523:   /* node 1 */
524:   s_u[2].i = i;
525:   s_u[2].j = j + 1;
526:   s_u[2].c = 0; /* Ux1 */
527:   s_u[3].i = i;
528:   s_u[3].j = j + 1;
529:   s_u[3].c = 1; /* Uy1 */

531:   /* node 2 */
532:   s_u[4].i = i + 1;
533:   s_u[4].j = j + 1;
534:   s_u[4].c = 0; /* Ux2 */
535:   s_u[5].i = i + 1;
536:   s_u[5].j = j + 1;
537:   s_u[5].c = 1; /* Uy2 */

539:   /* node 3 */
540:   s_u[6].i = i + 1;
541:   s_u[6].j = j;
542:   s_u[6].c = 0; /* Ux3 */
543:   s_u[7].i = i + 1;
544:   s_u[7].j = j;
545:   s_u[7].c = 1; /* Uy3 */
546:   return 0;
547: }

549: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
550: {
552:   /* get coords for the element */
553:   el_coords[NSD * 0 + 0] = _coords[ej][ei].x;
554:   el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
555:   el_coords[NSD * 1 + 0] = _coords[ej + 1][ei].x;
556:   el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
557:   el_coords[NSD * 2 + 0] = _coords[ej + 1][ei + 1].x;
558:   el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
559:   el_coords[NSD * 3 + 0] = _coords[ej][ei + 1].x;
560:   el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
561:   return 0;
562: }

564: static PetscErrorCode AssembleA_Elasticity(Mat A, DM elas_da, DM properties_da, Vec properties)
565: {
566:   DM                       cda;
567:   Vec                      coords;
568:   DMDACoor2d             **_coords;
569:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
570:   PetscInt                 sex, sey, mx, my;
571:   PetscInt                 ei, ej;
572:   PetscScalar              Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
573:   PetscScalar              el_coords[NODES_PER_EL * NSD];
574:   Vec                      local_properties;
575:   GaussPointCoefficients **props;
576:   PetscScalar             *prop_E, *prop_nu;

579:   /* setup for coords */
580:   DMGetCoordinateDM(elas_da, &cda);
581:   DMGetCoordinatesLocal(elas_da, &coords);
582:   DMDAVecGetArray(cda, coords, &_coords);

584:   /* setup for coefficients */
585:   DMCreateLocalVector(properties_da, &local_properties);
586:   DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties);
587:   DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties);
588:   DMDAVecGetArray(properties_da, local_properties, &props);

590:   DMDAGetElementsCorners(elas_da, &sex, &sey, 0);
591:   DMDAGetElementsSizes(elas_da, &mx, &my, 0);
592:   for (ej = sey; ej < sey + my; ej++) {
593:     for (ei = sex; ei < sex + mx; ei++) {
594:       /* get coords for the element */
595:       GetElementCoords(_coords, ei, ej, el_coords);

597:       /* get coefficients for the element */
598:       prop_E  = props[ej][ei].E;
599:       prop_nu = props[ej][ei].nu;

601:       /* initialise element stiffness matrix */
602:       PetscMemzero(Ae, sizeof(Ae));

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

607:       /* insert element matrix into global matrix */
608:       DMDAGetElementEqnums_u(u_eqn, ei, ej);
609:       MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
610:     }
611:   }
612:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
613:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

617:   DMDAVecRestoreArray(properties_da, local_properties, &props);
618:   VecDestroy(&local_properties);
619:   return 0;
620: }

622: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(ElasticityDOF **fields_F, MatStencil u_eqn[], PetscScalar Fe_u[])
623: {
624:   PetscInt n;

627:   for (n = 0; n < 4; n++) {
628:     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];
629:     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];
630:   }
631:   return 0;
632: }

634: static PetscErrorCode AssembleF_Elasticity(Vec F, DM elas_da, DM properties_da, Vec properties)
635: {
636:   DM                       cda;
637:   Vec                      coords;
638:   DMDACoor2d             **_coords;
639:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
640:   PetscInt                 sex, sey, mx, my;
641:   PetscInt                 ei, ej;
642:   PetscScalar              Fe[NODES_PER_EL * U_DOFS];
643:   PetscScalar              el_coords[NODES_PER_EL * NSD];
644:   Vec                      local_properties;
645:   GaussPointCoefficients **props;
646:   PetscScalar             *prop_fx, *prop_fy;
647:   Vec                      local_F;
648:   ElasticityDOF          **ff;

651:   /* setup for coords */
652:   DMGetCoordinateDM(elas_da, &cda);
653:   DMGetCoordinatesLocal(elas_da, &coords);
654:   DMDAVecGetArray(cda, coords, &_coords);

656:   /* setup for coefficients */
657:   DMGetLocalVector(properties_da, &local_properties);
658:   DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties);
659:   DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties);
660:   DMDAVecGetArray(properties_da, local_properties, &props);

662:   /* get access to the vector */
663:   DMGetLocalVector(elas_da, &local_F);
664:   VecZeroEntries(local_F);
665:   DMDAVecGetArray(elas_da, local_F, &ff);

667:   DMDAGetElementsCorners(elas_da, &sex, &sey, 0);
668:   DMDAGetElementsSizes(elas_da, &mx, &my, 0);
669:   for (ej = sey; ej < sey + my; ej++) {
670:     for (ei = sex; ei < sex + mx; ei++) {
671:       /* get coords for the element */
672:       GetElementCoords(_coords, ei, ej, el_coords);

674:       /* get coefficients for the element */
675:       prop_fx = props[ej][ei].fx;
676:       prop_fy = props[ej][ei].fy;

678:       /* initialise element stiffness matrix */
679:       PetscMemzero(Fe, sizeof(Fe));

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

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

687:       DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, Fe);
688:     }
689:   }

691:   DMDAVecRestoreArray(elas_da, local_F, &ff);
692:   DMLocalToGlobalBegin(elas_da, local_F, ADD_VALUES, F);
693:   DMLocalToGlobalEnd(elas_da, local_F, ADD_VALUES, F);
694:   DMRestoreLocalVector(elas_da, &local_F);

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

698:   DMDAVecRestoreArray(properties_da, local_properties, &props);
699:   DMRestoreLocalVector(properties_da, &local_properties);
700:   return 0;
701: }

703: static PetscErrorCode solve_elasticity_2d(PetscInt mx, PetscInt my)
704: {
705:   DM                       elas_da, da_prop;
706:   PetscInt                 u_dof, dof, stencil_width;
707:   Mat                      A;
708:   PetscInt                 mxl, myl;
709:   DM                       prop_cda, vel_cda;
710:   Vec                      prop_coords, vel_coords;
711:   PetscInt                 si, sj, nx, ny, i, j, p;
712:   Vec                      f, X;
713:   PetscInt                 prop_dof, prop_stencil_width;
714:   Vec                      properties, l_properties;
715:   MatNullSpace             matnull;
716:   PetscReal                dx, dy;
717:   PetscInt                 M, N;
718:   DMDACoor2d             **_prop_coords, **_vel_coords;
719:   GaussPointCoefficients **element_props;
720:   KSP                      ksp_E;
721:   PetscInt                 coefficient_structure = 0;
722:   PetscInt                 cpu_x, cpu_y, *lx = NULL, *ly = NULL;
723:   PetscBool                use_gp_coords = PETSC_FALSE;
724:   PetscBool                use_nonsymbc  = PETSC_FALSE;
725:   PetscBool                no_view       = PETSC_FALSE;
726:   PetscBool                flg;

729:   /* Generate the da for velocity and pressure */
730:   /*
731:    We use Q1 elements for the temperature.
732:    FEM has a 9-point stencil (BOX) or connectivity pattern
733:    Num nodes in each direction is mx+1, my+1
734:    */
735:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
736:   dof           = u_dof;
737:   stencil_width = 1;
738:   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);

740:   DMSetMatType(elas_da, MATAIJ);
741:   DMSetFromOptions(elas_da);
742:   DMSetUp(elas_da);

744:   DMDASetFieldName(elas_da, 0, "Ux");
745:   DMDASetFieldName(elas_da, 1, "Uy");

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

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

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

758:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
759:   prop_stencil_width = 0;
760:   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);
761:   DMSetFromOptions(da_prop);
762:   DMSetUp(da_prop);

764:   PetscFree(lx);
765:   PetscFree(ly);

767:   /* define centroid positions */
768:   DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
769:   dx = 1.0 / ((PetscReal)(M));
770:   dy = 1.0 / ((PetscReal)(N));

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

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

777:   DMCreateGlobalVector(da_prop, &properties);
778:   DMCreateLocalVector(da_prop, &l_properties);
779:   DMDAVecGetArray(da_prop, l_properties, &element_props);

781:   DMGetCoordinateDM(da_prop, &prop_cda);
782:   DMGetCoordinatesLocal(da_prop, &prop_coords);
783:   DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords);

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

787:   DMGetCoordinateDM(elas_da, &vel_cda);
788:   DMGetCoordinatesLocal(elas_da, &vel_coords);
789:   DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords);

791:   /* interpolate the coordinates */
792:   for (j = sj; j < sj + ny; j++) {
793:     for (i = si; i < si + nx; i++) {
794:       PetscInt    ngp;
795:       PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
796:       PetscScalar el_coords[8];

798:       GetElementCoords(_vel_coords, i, j, el_coords);
799:       ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

801:       for (p = 0; p < GAUSS_POINTS; p++) {
802:         PetscScalar gp_x, gp_y;
803:         PetscInt    n;
804:         PetscScalar xi_p[2], Ni_p[4];

806:         xi_p[0] = gp_xi[p][0];
807:         xi_p[1] = gp_xi[p][1];
808:         ConstructQ12D_Ni(xi_p, Ni_p);

810:         gp_x = 0.0;
811:         gp_y = 0.0;
812:         for (n = 0; n < NODES_PER_EL; n++) {
813:           gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
814:           gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
815:         }
816:         element_props[j][i].gp_coords[2 * p]     = gp_x;
817:         element_props[j][i].gp_coords[2 * p + 1] = gp_y;
818:       }
819:     }
820:   }

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

825:   for (j = sj; j < sj + ny; j++) {
826:     for (i = si; i < si + nx; i++) {
827:       PetscScalar              centroid_x = _prop_coords[j][i].x; /* centroids of cell */
828:       PetscScalar              centroid_y = _prop_coords[j][i].y;
829:       PETSC_UNUSED PetscScalar coord_x, coord_y;

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

834:         opts_E  = 1.0;
835:         opts_nu = 0.33;
836:         PetscOptionsGetScalar(NULL, NULL, "-iso_E", &opts_E, &flg);
837:         PetscOptionsGetScalar(NULL, NULL, "-iso_nu", &opts_nu, &flg);

839:         for (p = 0; p < GAUSS_POINTS; p++) {
840:           element_props[j][i].E[p]  = opts_E;
841:           element_props[j][i].nu[p] = opts_nu;

843:           element_props[j][i].fx[p] = 0.0;
844:           element_props[j][i].fy[p] = 0.0;
845:         }
846:       } else if (coefficient_structure == 1) { /* step */
847:         PetscScalar opts_E0, opts_nu0, opts_xc;
848:         PetscScalar opts_E1, opts_nu1;

850:         opts_E0 = opts_E1 = 1.0;
851:         opts_nu0 = opts_nu1 = 0.333;
852:         opts_xc             = 0.5;
853:         PetscOptionsGetScalar(NULL, NULL, "-step_E0", &opts_E0, &flg);
854:         PetscOptionsGetScalar(NULL, NULL, "-step_nu0", &opts_nu0, &flg);
855:         PetscOptionsGetScalar(NULL, NULL, "-step_E1", &opts_E1, &flg);
856:         PetscOptionsGetScalar(NULL, NULL, "-step_nu1", &opts_nu1, &flg);
857:         PetscOptionsGetScalar(NULL, NULL, "-step_xc", &opts_xc, &flg);

859:         for (p = 0; p < GAUSS_POINTS; p++) {
860:           coord_x = centroid_x;
861:           coord_y = centroid_y;
862:           if (use_gp_coords) {
863:             coord_x = element_props[j][i].gp_coords[2 * p];
864:             coord_y = element_props[j][i].gp_coords[2 * p + 1];
865:           }

867:           element_props[j][i].E[p]  = opts_E0;
868:           element_props[j][i].nu[p] = opts_nu0;
869:           if (PetscRealPart(coord_x) > PetscRealPart(opts_xc)) {
870:             element_props[j][i].E[p]  = opts_E1;
871:             element_props[j][i].nu[p] = opts_nu1;
872:           }

874:           element_props[j][i].fx[p] = 0.0;
875:           element_props[j][i].fy[p] = 0.0;
876:         }
877:       } else if (coefficient_structure == 2) { /* brick */
878:         PetscReal values_E[10];
879:         PetscReal values_nu[10];
880:         PetscInt  nbricks, maxnbricks;
881:         PetscInt  index, span;
882:         PetscInt  jj;

884:         flg        = PETSC_FALSE;
885:         maxnbricks = 10;
886:         PetscOptionsGetRealArray(NULL, NULL, "-brick_E", values_E, &maxnbricks, &flg);
887:         nbricks = maxnbricks;

890:         flg        = PETSC_FALSE;
891:         maxnbricks = 10;
892:         PetscOptionsGetRealArray(NULL, NULL, "-brick_nu", values_nu, &maxnbricks, &flg);

896:         span = 1;
897:         PetscOptionsGetInt(NULL, NULL, "-brick_span", &span, &flg);

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

904:         for (p = 0; p < GAUSS_POINTS; p++) {
905:           element_props[j][i].E[p]  = values_E[index];
906:           element_props[j][i].nu[p] = values_nu[index];
907:         }
908:       } else if (coefficient_structure == 3) { /* sponge */
909:         PetscScalar opts_E0, opts_nu0;
910:         PetscScalar opts_E1, opts_nu1;
911:         PetscInt    opts_t, opts_w;
912:         PetscInt    ii, jj, ci, cj;

914:         opts_E0 = opts_E1 = 1.0;
915:         opts_nu0 = opts_nu1 = 0.333;
916:         PetscOptionsGetScalar(NULL, NULL, "-sponge_E0", &opts_E0, &flg);
917:         PetscOptionsGetScalar(NULL, NULL, "-sponge_nu0", &opts_nu0, &flg);
918:         PetscOptionsGetScalar(NULL, NULL, "-sponge_E1", &opts_E1, &flg);
919:         PetscOptionsGetScalar(NULL, NULL, "-sponge_nu1", &opts_nu1, &flg);

921:         opts_t = opts_w = 1;
922:         PetscOptionsGetInt(NULL, NULL, "-sponge_t", &opts_t, &flg);
923:         PetscOptionsGetInt(NULL, NULL, "-sponge_w", &opts_w, &flg);

925:         ii = (i) / (opts_t + opts_w + opts_t);
926:         jj = (j) / (opts_t + opts_w + opts_t);

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

931:         for (p = 0; p < GAUSS_POINTS; p++) {
932:           element_props[j][i].E[p]  = opts_E0;
933:           element_props[j][i].nu[p] = opts_nu0;
934:         }
935:         if ((ci >= opts_t) && (ci < opts_t + opts_w)) {
936:           if ((cj >= opts_t) && (cj < opts_t + opts_w)) {
937:             for (p = 0; p < GAUSS_POINTS; p++) {
938:               element_props[j][i].E[p]  = opts_E1;
939:               element_props[j][i].nu[p] = opts_nu1;
940:             }
941:           }
942:         }
943:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
944:     }
945:   }
946:   DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords);

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

950:   DMDAVecRestoreArray(da_prop, l_properties, &element_props);
951:   DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties);
952:   DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties);

954:   PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL);
955:   if (!no_view) {
956:     DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for elasticity eqn.", "properties");
957:     DMDACoordViewGnuplot2d(elas_da, "mesh");
958:   }

960:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
961:   DMCreateMatrix(elas_da, &A);
962:   DMGetCoordinates(elas_da, &vel_coords);
963:   MatNullSpaceCreateRigidBody(vel_coords, &matnull);
964:   MatSetNearNullSpace(A, matnull);
965:   MatNullSpaceDestroy(&matnull);
966:   MatCreateVecs(A, &f, &X);

968:   /* assemble A11 */
969:   MatZeroEntries(A);
970:   VecZeroEntries(f);

972:   AssembleA_Elasticity(A, elas_da, da_prop, properties);
973:   /* build force vector */
974:   AssembleF_Elasticity(f, elas_da, da_prop, properties);

976:   KSPCreate(PETSC_COMM_WORLD, &ksp_E);
977:   KSPSetOptionsPrefix(ksp_E, "elas_"); /* elasticity */

979:   PetscOptionsGetBool(NULL, NULL, "-use_nonsymbc", &use_nonsymbc, &flg);
980:   /* solve */
981:   if (!use_nonsymbc) {
982:     Mat        AA;
983:     Vec        ff, XX;
984:     IS         is;
985:     VecScatter scat;

987:     DMDABCApplySymmetricCompression(elas_da, A, f, &is, &AA, &ff);
988:     VecDuplicate(ff, &XX);

990:     KSPSetOperators(ksp_E, AA, AA);
991:     KSPSetFromOptions(ksp_E);

993:     KSPSolve(ksp_E, ff, XX);

995:     /* push XX back into X */
996:     DMDABCApplyCompression(elas_da, NULL, X);

998:     VecScatterCreate(XX, NULL, X, is, &scat);
999:     VecScatterBegin(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD);
1000:     VecScatterEnd(scat, XX, X, INSERT_VALUES, SCATTER_FORWARD);
1001:     VecScatterDestroy(&scat);

1003:     MatDestroy(&AA);
1004:     VecDestroy(&ff);
1005:     VecDestroy(&XX);
1006:     ISDestroy(&is);
1007:   } else {
1008:     DMDABCApplyCompression(elas_da, A, f);

1010:     KSPSetOperators(ksp_E, A, A);
1011:     KSPSetFromOptions(ksp_E);

1013:     KSPSolve(ksp_E, f, X);
1014:   }

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

1019:   VecDestroy(&X);
1020:   VecDestroy(&f);
1021:   MatDestroy(&A);

1023:   DMDestroy(&elas_da);
1024:   DMDestroy(&da_prop);

1026:   VecDestroy(&properties);
1027:   VecDestroy(&l_properties);
1028:   return 0;
1029: }

1031: int main(int argc, char **args)
1032: {
1033:   PetscInt mx, my;

1036:   PetscInitialize(&argc, &args, (char *)0, help);
1037:   mx = my = 10;
1038:   PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL);
1039:   PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL);
1040:   solve_elasticity_2d(mx, my);
1041:   PetscFinalize();
1042:   return 0;
1043: }

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

1047: static PetscErrorCode BCApply_EAST(DM da, PetscInt d_idx, PetscScalar bc_val, Mat A, Vec b)
1048: {
1049:   DM                     cda;
1050:   Vec                    coords;
1051:   PetscInt               si, sj, nx, ny, i, j;
1052:   PetscInt               M, N;
1053:   DMDACoor2d           **_coords;
1054:   const PetscInt        *g_idx;
1055:   PetscInt              *bc_global_ids;
1056:   PetscScalar           *bc_vals;
1057:   PetscInt               nbcs;
1058:   PetscInt               n_dofs;
1059:   ISLocalToGlobalMapping ltogm;

1062:   /* enforce bc's */
1063:   DMGetLocalToGlobalMapping(da, &ltogm);
1064:   ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);

1066:   DMGetCoordinateDM(da, &cda);
1067:   DMGetCoordinatesLocal(da, &coords);
1068:   DMDAVecGetArray(cda, coords, &_coords);
1069:   DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
1070:   DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);

1072:   /* --- */

1074:   PetscMalloc1(ny * n_dofs, &bc_global_ids);
1075:   PetscMalloc1(ny * n_dofs, &bc_vals);

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

1080:   i = nx - 1;
1081:   for (j = 0; j < ny; j++) {
1082:     PetscInt                 local_id;
1083:     PETSC_UNUSED PetscScalar coordx, coordy;

1085:     local_id = i + j * nx;

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

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

1092:     bc_vals[j] = bc_val;
1093:   }
1094:   ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1095:   nbcs = 0;
1096:   if ((si + nx) == (M)) nbcs = ny;

1098:   if (b) {
1099:     VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1100:     VecAssemblyBegin(b);
1101:     VecAssemblyEnd(b);
1102:   }
1103:   if (A) MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0);

1105:   PetscFree(bc_vals);
1106:   PetscFree(bc_global_ids);

1108:   DMDAVecRestoreArray(cda, coords, &_coords);
1109:   return 0;
1110: }

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

1127:   /* enforce bc's */
1128:   DMGetLocalToGlobalMapping(da, &ltogm);
1129:   ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);

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

1137:   /* --- */

1139:   PetscMalloc1(ny * n_dofs, &bc_global_ids);
1140:   PetscMalloc1(ny * n_dofs, &bc_vals);

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

1145:   i = 0;
1146:   for (j = 0; j < ny; j++) {
1147:     PetscInt                 local_id;
1148:     PETSC_UNUSED PetscScalar coordx, coordy;

1150:     local_id = i + j * nx;

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

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

1157:     bc_vals[j] = bc_val;
1158:   }
1159:   ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1160:   nbcs = 0;
1161:   if (si == 0) nbcs = ny;

1163:   if (b) {
1164:     VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1165:     VecAssemblyBegin(b);
1166:     VecAssemblyEnd(b);
1167:   }
1168:   if (A) MatZeroRows(A, nbcs, bc_global_ids, 1.0, 0, 0);

1170:   PetscFree(bc_vals);
1171:   PetscFree(bc_global_ids);

1173:   DMDAVecRestoreArray(cda, coords, &_coords);
1174:   return 0;
1175: }

1177: static PetscErrorCode DMDABCApplyCompression(DM elas_da, Mat A, Vec f)
1178: {
1180:   BCApply_EAST(elas_da, 0, -1.0, A, f);
1181:   BCApply_EAST(elas_da, 1, 0.0, A, f);
1182:   BCApply_WEST(elas_da, 0, 1.0, A, f);
1183:   BCApply_WEST(elas_da, 1, 0.0, A, f);
1184:   return 0;
1185: }

1187: static PetscErrorCode Orthogonalize(PetscInt n, Vec *vecs)
1188: {
1189:   PetscInt    i, j;
1190:   PetscScalar dot;

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

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

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

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

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

1240:   ISCreateGeneral(PETSC_COMM_WORLD, cnt, unconstrained, PETSC_COPY_VALUES, &is);
1241:   PetscFree(unconstrained);
1242:   ISSetBlockSize(is, 2);

1244:   /* define correction for dirichlet in the rhs */
1245:   MatMult(A, x, f);
1246:   VecScale(f, -1.0);

1248:   /* get new matrix */
1249:   MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, AA);
1250:   /* get new vector */
1251:   MatCreateVecs(*AA, NULL, ff);

1253:   VecScatterCreate(f, is, *ff, NULL, &scat);
1254:   VecScatterBegin(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD);
1255:   VecScatterEnd(scat, f, *ff, INSERT_VALUES, SCATTER_FORWARD);

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

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

1278:   VecScatterDestroy(&scat);

1280:   *dofs = is;
1281:   VecDestroy(&x);
1282:   return 0;
1283: }

1285: /*TEST

1287:    build:
1288:       requires: !complex !single

1290:    test:
1291:       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
1292:       output_file: output/ex49_1.out

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

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

1304:    test:
1305:       suffix: 4
1306:       nsize: 4
1307:       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
1308:       requires: ml

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

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

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

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

1330:    test:
1331:       suffix: hypre_nullspace
1332:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
1333:       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

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

1340:    test:
1341:       nsize: 4
1342:       suffix: bddc_unsym
1343:       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

1345:    test:
1346:       nsize: 4
1347:       suffix: bddc_unsym_deluxe
1348:       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

1350:    test:
1351:       nsize: 4
1352:       suffix: fetidp_unsym_deluxe
1353:       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

1355:    test:
1356:       nsize: 4
1357:       suffix: bddc_layerjump
1358:       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

1360:    test:
1361:       nsize: 4
1362:       suffix: bddc_subdomainjump
1363:       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

1365:    test:
1366:       nsize: 9
1367:       suffix: bddc_subdomainjump_deluxe
1368:       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
1369: TEST*/