Actual source code: ex56.c

petsc-master 2019-09-15
Report Typos and Errors
  1: static char help[] = "3D, tri-quadratic hexahedra (Q1), displacement finite element formulation\n\
  2: of linear elasticity.  E=1.0, nu=1/3.\n\
  3: Unit cube domain with Dirichlet boundary\n\n";

  5:  #include <petscdmplex.h>
  6:  #include <petscsnes.h>
  7:  #include <petscds.h>
  8:  #include <petscdmforest.h>

 10: static PetscReal s_soft_alpha=1.e-3;
 11: static PetscReal s_mu=0.4;
 12: static PetscReal s_lambda=0.4;

 14: static void f0_bd_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 15:                        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 16:                        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 17:                        PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 18: {
 19:   f0[0] = 1;     /* x direction pull */
 20:   f0[1] = -x[2]; /* add a twist around x-axis */
 21:   f0[2] =  x[1];
 22: }

 24: static void f1_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 25:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 26:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 27:                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 28: {
 29:   const PetscInt Ncomp = dim;
 30:   PetscInt       comp, d;
 31:   for (comp = 0; comp < Ncomp; ++comp) {
 32:     for (d = 0; d < dim; ++d) {
 33:       f1[comp*dim+d] = 0.0;
 34:     }
 35:   }
 36: }

 38: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
 39: static void f1_u_3d_alpha(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 40:                           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 41:                           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 42:                           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 43: {
 44:   PetscReal trace,mu=s_mu,lambda=s_lambda,rad;
 45:   PetscInt i,j;
 46:   for (i=0,rad=0.;i<dim;i++) {
 47:     PetscReal t=x[i];
 48:     rad += t*t;
 49:   }
 50:   rad = PetscSqrtReal(rad);
 51:   if (rad>0.25) {
 52:     mu *= s_soft_alpha;
 53:     lambda *= s_soft_alpha; /* we could keep the bulk the same like rubberish */
 54:   }
 55:   for (i=0,trace=0; i < dim; ++i) {
 56:     trace += PetscRealPart(u_x[i*dim+i]);
 57:   }
 58:   for (i=0; i < dim; ++i) {
 59:     for (j=0; j < dim; ++j) {
 60:       f1[i*dim+j] = mu*(u_x[i*dim+j]+u_x[j*dim+i]);
 61:     }
 62:     f1[i*dim+i] += lambda*trace;
 63:   }
 64: }

 66: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
 67: static void f1_u_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 68:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 69:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 70:                     PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 71: {
 72:   PetscReal trace,mu=s_mu,lambda=s_lambda;
 73:   PetscInt i,j;
 74:   for (i=0,trace=0; i < dim; ++i) {
 75:     trace += PetscRealPart(u_x[i*dim+i]);
 76:   }
 77:   for (i=0; i < dim; ++i) {
 78:     for (j=0; j < dim; ++j) {
 79:       f1[i*dim+j] = mu*(u_x[i*dim+j]+u_x[j*dim+i]);
 80:     }
 81:     f1[i*dim+i] += lambda*trace;
 82:   }
 83: }

 85: /* 3D elasticity */
 86: #define IDX(ii,jj,kk,ll) (27*ii+9*jj+3*kk+ll)

 88: void g3_uu_3d_private( PetscScalar g3[], const PetscReal mu, const PetscReal lambda)
 89: {
 90:   if (1) {
 91:     g3[0] += lambda;
 92:     g3[0] += mu;
 93:     g3[0] += mu;
 94:     g3[4] += lambda;
 95:     g3[8] += lambda;
 96:     g3[10] += mu;
 97:     g3[12] += mu;
 98:     g3[20] += mu;
 99:     g3[24] += mu;
100:     g3[28] += mu;
101:     g3[30] += mu;
102:     g3[36] += lambda;
103:     g3[40] += lambda;
104:     g3[40] += mu;
105:     g3[40] += mu;
106:     g3[44] += lambda;
107:     g3[50] += mu;
108:     g3[52] += mu;
109:     g3[56] += mu;
110:     g3[60] += mu;
111:     g3[68] += mu;
112:     g3[70] += mu;
113:     g3[72] += lambda;
114:     g3[76] += lambda;
115:     g3[80] += lambda;
116:     g3[80] += mu;
117:     g3[80] += mu;
118:   } else {
119:     int        i,j,k,l;
120:     static int cc=-1;
121:     cc++;
122:     for (i = 0; i < 3; ++i) {
123:       for (j = 0; j < 3; ++j) {
124:         for (k = 0; k < 3; ++k) {
125:           for (l = 0; l < 3; ++l) {
126:             if (k==l && i==j) g3[IDX(i,j,k,l)] += lambda;
127:             if (i==k && j==l) g3[IDX(i,j,k,l)] += mu;
128:             if (i==l && j==k) g3[IDX(i,j,k,l)] += mu;
129:             if (k==l && i==j && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += lambda;\n",IDX(i,j,k,l));
130:             if (i==k && j==l && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += mu;\n",IDX(i,j,k,l));
131:             if (i==l && j==k && !cc) (void) PetscPrintf(PETSC_COMM_WORLD,"g3[%d] += mu;\n",IDX(i,j,k,l));
132:           }
133:         }
134:       }
135:     }
136:   }
137: }

139: static void g3_uu_3d_alpha(PetscInt dim, PetscInt Nf, PetscInt NfAux,
140:                            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
141:                            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
142:                            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
143: {
144:   PetscReal mu=s_mu, lambda=s_lambda,rad;
145:   PetscInt i;
146:   for (i=0,rad=0.;i<dim;i++) {
147:     PetscReal t=x[i];
148:     rad += t*t;
149:   }
150:   rad = PetscSqrtReal(rad);
151:   if (rad>0.25) {
152:     mu *= s_soft_alpha;
153:     lambda *= s_soft_alpha; /* we could keep the bulk the same like rubberish */
154:   }
155:   g3_uu_3d_private(g3,mu,lambda);
156: }

158: static void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
159:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
160:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
161:                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
162: {
163:   g3_uu_3d_private(g3,s_mu,s_lambda);
164: }

166: static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
167:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
168:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
169:                  PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
170: {
171:   const    PetscInt Ncomp = dim;
172:   PetscInt comp;

174:   for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 0.0;
175: }

177: /* PI_i (x_i^4 - x_i^2) */
178: static void f0_u_x4(PetscInt dim, PetscInt Nf, PetscInt NfAux,
179:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
180:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
181:                     PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
182: {
183:   const    PetscInt Ncomp = dim;
184:   PetscInt comp,i;

186:   for (comp = 0; comp < Ncomp; ++comp) {
187:     f0[comp] = 1e5;
188:     for (i = 0; i < Ncomp; ++i) {
189:       f0[comp] *= /* (comp+1)* */(x[i]*x[i]*x[i]*x[i] - x[i]*x[i]); /* assumes (0,1]^D domain */
190:     }
191:   }
192: }

194: PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
195: {
196:   const PetscInt Ncomp = dim;
197:   PetscInt       comp;

199:   for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0;
200:   return 0;
201: }

203: int main(int argc,char **args)
204: {
205:   Mat            Amat;
207:   SNES           snes;
208:   KSP            ksp;
209:   MPI_Comm       comm;
210:   PetscMPIInt    rank;
211:   PetscLogStage  stage[7];
212:   PetscBool      test_nonzero_cols=PETSC_FALSE,use_nearnullspace=PETSC_TRUE,attach_nearnullspace=PETSC_FALSE;
213:   Vec            xx,bb;
214:   PetscInt       iter,i,N,dim=3,cells[3]={1,1,1},max_conv_its,local_sizes[7],run_type=1;
215:   DM             dm,distdm,basedm;
216:   PetscBool      flg;
217:   char           convType[256];
218:   PetscReal      Lx,mdisp[10],err[10];
219:   const char * const options[10] = {"-ex56_dm_refine 0",
220:                                     "-ex56_dm_refine 1",
221:                                     "-ex56_dm_refine 2",
222:                                     "-ex56_dm_refine 3",
223:                                     "-ex56_dm_refine 4",
224:                                     "-ex56_dm_refine 5",
225:                                     "-ex56_dm_refine 6",
226:                                     "-ex56_dm_refine 7",
227:                                     "-ex56_dm_refine 8",
228:                                     "-ex56_dm_refine 9"};
230:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
231:   comm = PETSC_COMM_WORLD;
232:   MPI_Comm_rank(comm, &rank);
233:   /* options */
234:   PetscOptionsBegin(comm,NULL,"3D bilinear Q1 elasticity options","");
235:   {
236:     i = 3;
237:     PetscOptionsIntArray("-cells", "Number of (flux tube) processor in each dimension", "ex56.c", cells, &i, NULL);

239:     Lx = 1.; /* or ne for rod */
240:     max_conv_its = 3;
241:     PetscOptionsInt("-max_conv_its","Number of iterations in convergence study","",max_conv_its,&max_conv_its,NULL);
242:     if (max_conv_its<=0 || max_conv_its>7) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_USER, "Bad number of iterations for convergence test (%D)",max_conv_its);
243:     PetscOptionsReal("-lx","Length of domain","",Lx,&Lx,NULL);
244:     PetscOptionsReal("-alpha","material coefficient inside circle","",s_soft_alpha,&s_soft_alpha,NULL);
245:     PetscOptionsBool("-test_nonzero_cols","nonzero test","",test_nonzero_cols,&test_nonzero_cols,NULL);
246:     PetscOptionsBool("-use_mat_nearnullspace","MatNearNullSpace API test","",use_nearnullspace,&use_nearnullspace,NULL);
247:     PetscOptionsBool("-attach_mat_nearnullspace","MatNearNullSpace API test (via MatSetNearNullSpace)","",attach_nearnullspace,&attach_nearnullspace,NULL);
248:     PetscOptionsInt("-run_type","0: twisting load on cantalever, 1: 3rd order accurate convergence test","",run_type,&run_type,NULL);
249:     i = 3;
250:     PetscOptionsInt("-mat_block_size","","",i,&i,&flg);
251:     if (!flg || i!=3) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_USER, "'-mat_block_size 3' must be set (%D) and = 3 (%D)",flg,flg? i : 3);
252:   }
253:   PetscOptionsEnd();
254:   PetscLogStageRegister("Mesh Setup", &stage[6]);
255:   PetscLogStageRegister("1st Setup", &stage[0]);
256:   PetscLogStageRegister("1st Solve", &stage[1]);

258:   /* create DM, Plex calls DMSetup */
259:   PetscLogStagePush(stage[6]);
260:   DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, cells, NULL, NULL, NULL, PETSC_TRUE, &dm);
261:   {
262:     DMLabel         label;
263:     IS              is;
264:     DMCreateLabel(dm, "boundary");
265:     DMGetLabel(dm, "boundary", &label);
266:     DMPlexMarkBoundaryFaces(dm, 1, label);
267:     if (run_type==0) {
268:       DMGetStratumIS(dm, "boundary", 1,  &is);
269:       DMCreateLabel(dm,"Faces");
270:       if (is) {
271:         PetscInt        d, f, Nf;
272:         const PetscInt *faces;
273:         PetscInt        csize;
274:         PetscSection    cs;
275:         Vec             coordinates ;
276:         DM              cdm;
277:         ISGetLocalSize(is, &Nf);
278:         ISGetIndices(is, &faces);
279:         DMGetCoordinatesLocal(dm, &coordinates);
280:         DMGetCoordinateDM(dm, &cdm);
281:         DMGetLocalSection(cdm, &cs);
282:         /* Check for each boundary face if any component of its centroid is either 0.0 or 1.0 */
283:         for (f = 0; f < Nf; ++f) {
284:           PetscReal   faceCoord;
285:           PetscInt    b,v;
286:           PetscScalar *coords = NULL;
287:           PetscInt    Nv;
288:           DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
289:           Nv   = csize/dim; /* Calculate mean coordinate vector */
290:           for (d = 0; d < dim; ++d) {
291:             faceCoord = 0.0;
292:             for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v*dim+d]);
293:             faceCoord /= Nv;
294:             for (b = 0; b < 2; ++b) {
295:               if (PetscAbs(faceCoord - b) < PETSC_SMALL) { /* domain have not been set yet, still [0,1]^3 */
296:                 DMSetLabelValue(dm, "Faces", faces[f], d*2+b+1);
297:               }
298:             }
299:           }
300:           DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
301:         }
302:         ISRestoreIndices(is, &faces);
303:       }
304:       ISDestroy(&is);
305:       DMGetLabel(dm, "Faces", &label);
306:       DMPlexLabelComplete(dm, label);
307:     }
308:   }
309:   {
310:     PetscInt    dimEmbed, i;
311:     PetscInt    nCoords;
312:     PetscScalar *coords,bounds[] = {0,1,-.5,.5,-.5,.5,}; /* x_min,x_max,y_min,y_max */
313:     Vec         coordinates;
314:     bounds[1] = Lx;
315:     if (run_type==1) {
316:       for (i = 0; i < 2*dim; i++) bounds[i] = (i%2) ? 1 : 0;
317:     }
318:     DMGetCoordinatesLocal(dm,&coordinates);
319:     DMGetCoordinateDim(dm,&dimEmbed);
320:     if (dimEmbed != dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"dimEmbed != dim %D",dimEmbed);
321:     VecGetLocalSize(coordinates,&nCoords);
322:     if (nCoords % dimEmbed) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Coordinate vector the wrong size");
323:     VecGetArray(coordinates,&coords);
324:     for (i = 0; i < nCoords; i += dimEmbed) {
325:       PetscInt    j;
326:       PetscScalar *coord = &coords[i];
327:       for (j = 0; j < dimEmbed; j++) {
328:         coord[j] = bounds[2 * j] + coord[j] * (bounds[2 * j + 1] - bounds[2 * j]);
329:       }
330:     }
331:     VecRestoreArray(coordinates,&coords);
332:     DMSetCoordinatesLocal(dm,coordinates);
333:   }

335:   /* convert to p4est, and distribute */

337:   PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
338:   PetscOptionsFList("-dm_type","Convert DMPlex to another format (should not be Plex!)","ex56.c",DMList,DMPLEX,convType,256,&flg);
339:   PetscOptionsEnd();
340:   if (flg) {
341:     DM newdm;
342:     DMConvert(dm,convType,&newdm);
343:     if (newdm) {
344:       const char *prefix;
345:       PetscBool isForest;
346:       PetscObjectGetOptionsPrefix((PetscObject)dm,&prefix);
347:       PetscObjectSetOptionsPrefix((PetscObject)newdm,prefix);
348:       DMIsForest(newdm,&isForest);
349:       if (isForest) {
350:       } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Converted to non Forest?");
351:       DMDestroy(&dm);
352:       dm   = newdm;
353:     } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Convert failed?");
354:   } else {
355:     PetscPartitioner part;
356:     /* Plex Distribute mesh over processes */
357: #if 1
358:     DMPlexGetPartitioner(dm,&part);
359:     PetscPartitionerSetFromOptions(part);
360: #endif
361:     DMPlexDistribute(dm, 0, NULL, &distdm);
362:     if (distdm) {
363:       const char *prefix;
364:       PetscObjectGetOptionsPrefix((PetscObject)dm,&prefix);
365:       PetscObjectSetOptionsPrefix((PetscObject)distdm,prefix);
366:       DMDestroy(&dm);
367:       dm   = distdm;
368:     }
369:   }
370:   PetscLogStagePop();
371:   basedm = dm; dm = NULL;

373:   for (iter=0 ; iter<max_conv_its ; iter++) {
374:     PetscLogStagePush(stage[6]);
375:     /* make new DM */
376:     DMClone(basedm, &dm);
377:     PetscObjectSetOptionsPrefix((PetscObject) dm, "ex56_");
378:     PetscObjectSetName( (PetscObject)dm,"Mesh");
379:     PetscOptionsClearValue(NULL,"-ex56_dm_refine");
380:     PetscOptionsInsertString(NULL,options[iter]);
381:     DMSetFromOptions(dm); /* refinement done here in Plex, p4est */
382:     /* snes */
383:     SNESCreate(comm, &snes);
384:     SNESSetDM(snes, dm);
385:     /* fem */
386:     {
387:       const PetscInt Ncomp = dim;
388:       const PetscInt components[] = {0,1,2};
389:       const PetscInt Nfid = 1, Npid = 1;
390:       const PetscInt fid[] = {1}; /* The fixed faces (x=0) */
391:       const PetscInt pid[] = {2}; /* The faces with loading (x=L_x) */
392:       PetscFE        fe;
393:       PetscDS        prob;
394:       DM             cdm = dm;

396:       PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, dim, PETSC_FALSE, NULL, PETSC_DECIDE, &fe); /* elasticity */
397:       PetscObjectSetName((PetscObject) fe, "deformation");
398:       /* FEM prob */
399:       DMSetField(dm, 0, NULL, (PetscObject) fe);
400:       DMCreateDS(dm);
401:       DMGetDS(dm, &prob);
402:       /* setup problem */
403:       if (run_type==1) {
404:         PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d);
405:         PetscDSSetResidual(prob, 0, f0_u_x4, f1_u_3d);
406:       } else {
407:         PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu_3d_alpha);
408:         PetscDSSetResidual(prob, 0, f0_u, f1_u_3d_alpha);
409:         PetscDSSetBdResidual(prob, 0, f0_bd_u_3d, f1_bd_u);
410:       }
411:       /* bcs */
412:       if (run_type==1) {
413:         PetscInt id = 1;
414:         DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "boundary", 0, 0, NULL, (void (*)(void)) zero, 1, &id, NULL);
415:       } else {
416:         PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "fixed", "Faces", 0, Ncomp, components, (void (*)(void)) zero, Nfid, fid, NULL);
417:         PetscDSAddBoundary(prob, DM_BC_NATURAL, "traction", "Faces", 0, Ncomp, components, NULL, Npid, pid, NULL);
418:       }
419:       while (cdm) {
420:         DMCopyDisc(dm, cdm);
421:         DMGetCoarseDM(cdm, &cdm);
422:       }
423:       PetscFEDestroy(&fe);
424:     }
425:     /* vecs & mat */
426:     DMCreateGlobalVector(dm,&xx);
427:     VecDuplicate(xx, &bb);
428:     PetscObjectSetName((PetscObject) bb, "b");
429:     PetscObjectSetName((PetscObject) xx, "u");
430:     DMCreateMatrix(dm, &Amat);
431:     VecGetSize(bb,&N);
432:     local_sizes[iter] = N;
433:     PetscPrintf(PETSC_COMM_WORLD,"[%d] %D global equations, %D vertices\n",rank,N,N/dim);
434:     if ((use_nearnullspace || attach_nearnullspace) && N/dim > 1) {
435:       /* Set up the near null space (a.k.a. rigid body modes) that will be used by the multigrid preconditioner */
436:       DM           subdm;
437:       MatNullSpace nearNullSpace;
438:       PetscInt     fields = 0;
439:       PetscObject  deformation;
440:       DMCreateSubDM(dm, 1, &fields, NULL, &subdm);
441:       DMPlexCreateRigidBody(subdm, &nearNullSpace);
442:       DMGetField(dm, 0, NULL, &deformation);
443:       PetscObjectCompose(deformation, "nearnullspace", (PetscObject) nearNullSpace);
444:       DMDestroy(&subdm);
445:       if (attach_nearnullspace) {
446:         MatSetNearNullSpace(Amat,nearNullSpace);
447:       }
448:       MatNullSpaceDestroy(&nearNullSpace); /* created by DM and destroyed by Mat */
449:     }
450:     DMPlexSetSNESLocalFEM(dm,NULL,NULL,NULL);
451:     SNESSetJacobian(snes, Amat, Amat, NULL, NULL);
452:     SNESSetFromOptions(snes);
453:     DMSetUp(dm);
454:     PetscLogStagePop();
455:     PetscLogStagePush(stage[0]);
456:     /* ksp */
457:     SNESGetKSP(snes, &ksp);
458:     KSPSetComputeSingularValues(ksp,PETSC_TRUE);
459:     /* test BCs */
460:     VecZeroEntries(xx);
461:     if (test_nonzero_cols) {
462:       if (rank==0) VecSetValue(xx,0,1.0,INSERT_VALUES);
463:       VecAssemblyBegin(xx);
464:       VecAssemblyEnd(xx);
465:     }
466:     VecZeroEntries(bb);
467:     VecGetSize(bb,&i);
468:     local_sizes[iter] = i;
469:     PetscPrintf(PETSC_COMM_WORLD,"[%d] %D equations in vector, %D vertices\n",rank,i,i/dim);
470:     /* setup solver, dummy solve to really setup */
471:     if (0) {
472:       KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);
473:       SNESSolve(snes, bb, xx);
474:       KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,50);
475:       VecZeroEntries(xx);
476:     }
477:     PetscLogStagePop();
478:     /* solve */
479:     PetscLogStagePush(stage[1]);
480:     SNESSolve(snes, bb, xx);
481:     PetscLogStagePop();
482:     VecNorm(xx,NORM_INFINITY,&mdisp[iter]);
483:     DMViewFromOptions(dm, NULL, "-dm_view");
484:     {
485:       PetscViewer       viewer = NULL;
486:       PetscViewerFormat fmt;
487:       PetscOptionsGetViewer(comm,NULL,"ex56_","-vec_view",&viewer,&fmt,&flg);
488:       if (flg) {
489:         PetscViewerPushFormat(viewer,fmt);
490:         VecView(xx,viewer);
491:         VecView(bb,viewer);
492:         PetscViewerPopFormat(viewer);
493:       }
494:       PetscViewerDestroy(&viewer);
495:     }
496:     /* Free work space */
497:     DMDestroy(&dm);
498:     SNESDestroy(&snes);
499:     VecDestroy(&xx);
500:     VecDestroy(&bb);
501:     MatDestroy(&Amat);
502:   }
503:   DMDestroy(&basedm);
504:   if (run_type==1) {
505:     err[0] = 59.975208 - mdisp[0]; /* error with what I think is the exact solution */
506:   } else {
507:     err[0] = 171.038 - mdisp[0];
508:   }
509:   for (iter=1 ; iter<max_conv_its ; iter++) {
510:     if (run_type==1) {
511:       err[iter] = 59.975208 - mdisp[iter];
512:     } else {
513:       err[iter] = 171.038 - mdisp[iter];
514:     }
515:     PetscPrintf(PETSC_COMM_WORLD,"[%d] %D) N=%12D, max displ=%9.7e, disp diff=%9.2e, error=%4.3e, rate=%3.2g\n",rank,iter,local_sizes[iter],mdisp[iter],
516:                 mdisp[iter]-mdisp[iter-1],err[iter],PetscLogReal(err[iter-1]/err[iter])/PetscLogReal(2.));
517:   }

519:   PetscFinalize();
520:   return ierr;
521: }

523: /*TEST

525:   test:
526:     suffix: 0
527:     nsize: 4
528:     requires: !single
529:     args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 1000 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -ksp_converged_reason -snes_monitor_short -ksp_monitor_short -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi -petscpartitioner_type simple -mat_block_size 3 -matptap_via scalable -ex56_dm_view -run_type 1
530:     timeoutfactor: 2

532:   # HYPRE PtAP broken with complex numbers
533:   test:
534:     suffix: hypre
535:     requires: hypre !single !complex
536:     nsize: 4
537:     args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -ex56_dm_refine 1 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -ksp_converged_reason -use_mat_nearnullspace true -mat_block_size 3 -petscpartitioner_type simple

539:   test:
540:     suffix: ml
541:     requires: ml !single
542:     nsize: 4
543:     args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -ex56_dm_refine 1 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_converged_reason -ksp_rtol 1.e-8 -pc_type ml -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 3 -mg_levels_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type sor -mg_levels_esteig_ksp_type cg -mat_block_size 3 -petscpartitioner_type simple -use_mat_nearnullspace

545:   test:
546:     nsize: 4
547:     requires: parmetis !single
548:     args: -cells 2,2,1 -max_conv_its 2 -petscspace_degree 2 -snes_max_it 2 -ksp_max_it 100 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -snes_rtol 1.e-10 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -snes_converged_reason -use_mat_nearnullspace true -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mg_levels_pc_type jacobi -pc_gamg_mat_partitioning_type parmetis -mat_block_size 3 -run_type 1 -pc_gamg_repartition true

550:   test:
551:     suffix: bddc
552:     nsize: 4
553:     requires: !single
554:     args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -ex56_dm_refine 1 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -mat_block_size 3 -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type {{sbaij baij aij}} -pc_type bddc

556:   testset:
557:     nsize: 4
558:     requires: !single
559:     args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -ex56_dm_refine 1 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-10 -ksp_converged_reason -mat_block_size 3 -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type aij -pc_type bddc -attach_mat_nearnullspace {{0 1}separate output}
560:     test:
561:       suffix: bddc_approx_gamg
562:       args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -approximate -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -prefix_pop -prefix_push pc_bddc_neumann_ -approximate -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -pc_gamg_coarse_eq_limit 10 -pc_gamg_reuse_interpolation true -pc_gamg_square_graph 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale .0 -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -prefix_pop
563:     # HYPRE PtAP broken with complex numbers
564:     test:
565:       requires: hypre !complex
566:       suffix: bddc_approx_hypre
567:       args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -pc_type hypre -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_strong_threshold 0.75 -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -prefix_pop -prefix_push pc_bddc_neumann_ -pc_type hypre -pc_hypre_boomeramg_no_CF true -pc_hypre_boomeramg_strong_threshold 0.75 -pc_hypre_boomeramg_agg_nl 1 -pc_hypre_boomeramg_coarsen_type HMIS -pc_hypre_boomeramg_interp_type ext+i -prefix_pop
568:     test:
569:       requires: ml
570:       suffix: bddc_approx_ml
571:       args: -pc_bddc_switch_static -prefix_push pc_bddc_dirichlet_ -approximate -pc_type ml -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -prefix_pop -prefix_push pc_bddc_neumann_ -approximate -pc_type ml -mg_levels_ksp_max_it 1 -mg_levels_ksp_type chebyshev -mg_levels_esteig_ksp_type cg -mg_levels_esteig_ksp_max_it 10 -prefix_pop

573:   test:
574:     suffix: fetidp
575:     nsize: 4
576:     requires: !single
577:     args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -ex56_dm_refine 1 -petscspace_degree 2 -ksp_type fetidp -fetidp_ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -mat_block_size 3 -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type {{sbaij baij aij}}

579:   test:
580:     suffix: bddc_elast
581:     nsize: 4
582:     requires: !single
583:     args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -ex56_dm_refine 1 -petscspace_degree 2 -ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -mat_block_size 3 -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type sbaij -pc_type bddc -pc_bddc_monolithic -attach_mat_nearnullspace

585:   test:
586:     suffix: fetidp_elast
587:     nsize: 4
588:     requires: !single
589:     args: -cells 2,2,1 -max_conv_its 2 -lx 1. -alpha .01 -ex56_dm_refine 1 -petscspace_degree 2 -ksp_type fetidp -fetidp_ksp_type cg -ksp_monitor_short -ksp_rtol 1.e-8 -ksp_converged_reason -mat_block_size 3 -petscpartitioner_type simple -ex56_dm_mat_type is -matis_localmat_type sbaij -fetidp_bddc_pc_bddc_monolithic -attach_mat_nearnullspace

591: TEST*/