Actual source code: dmplexsnes.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/snesimpl.h>
  3: #include <petscds.h>
  4: #include <petsc/private/petscimpl.h>
  5: #include <petsc/private/petscfeimpl.h>

  7: static void pressure_Private(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar p[])
  8: {
  9:   p[0] = u[uOff[1]];
 10: }

 12: /*
 13:   SNESCorrectDiscretePressure_Private - Add a vector in the nullspace to make the continuum integral of the pressure field equal to zero. This is normally used only to evaluate convergence rates for the pressure accurately.

 15:   Collective

 17:   Input Parameters:
 18: + snes      - The SNES
 19: . pfield    - The field number for pressure
 20: . nullspace - The pressure nullspace
 21: . u         - The solution vector
 22: - ctx       - An optional user context

 24:   Output Parameter:
 25: . u         - The solution with a continuum pressure integral of zero

 27:   Notes:
 28:   If int(u) = a and int(n) = b, then int(u - a/b n) = a - a/b b = 0. We assume that the nullspace is a single vector given explicitly.

 30:   Level: developer

 32: .seealso: `SNESConvergedCorrectPressure()`
 33: */
 34: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *ctx)
 35: {
 36:   DM          dm;
 37:   PetscDS     ds;
 38:   const Vec  *nullvecs;
 39:   PetscScalar pintd, *intc, *intn;
 40:   MPI_Comm    comm;
 41:   PetscInt    Nf, Nv;

 43:   PetscObjectGetComm((PetscObject)snes, &comm);
 44:   SNESGetDM(snes, &dm);
 47:   DMGetDS(dm, &ds);
 48:   PetscDSSetObjective(ds, pfield, pressure_Private);
 49:   MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs);
 51:   VecDot(nullvecs[0], u, &pintd);
 53:   PetscDSGetNumFields(ds, &Nf);
 54:   PetscMalloc2(Nf, &intc, Nf, &intn);
 55:   DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx);
 56:   DMPlexComputeIntegralFEM(dm, u, intc, ctx);
 57:   VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0]);
 58: #if defined(PETSC_USE_DEBUG)
 59:   DMPlexComputeIntegralFEM(dm, u, intc, ctx);
 61: #endif
 62:   PetscFree2(intc, intn);
 63:   return 0;
 64: }

 66: /*@C
 67:    SNESConvergedCorrectPressure - Convergence test that adds a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
 68:    This is normally used only to evaluate convergence rates for the pressure accurately. The convergence test itself just mimics `SNESConvergedDefault()`.

 70:    Logically Collective

 72:    Input Parameters:
 73: +  snes - the `SNES` context
 74: .  it - the iteration (0 indicates before any Newton steps)
 75: .  xnorm - 2-norm of current iterate
 76: .  snorm - 2-norm of current step
 77: .  fnorm - 2-norm of function at current iterate
 78: -  ctx   - Optional user context

 80:    Output Parameter:
 81: .  reason  - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FNORM_NAN`

 83:    Notes:
 84:    In order to use this convergence test, you must set up several PETSc structures. First fields must be added to the `DM`, and a `PetscDS`
 85:    must be created with discretizations of those fields. We currently assume that the pressure field has index 1.
 86:    The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface.
 87:    Last we must be able to integrate the pressure over the domain, so the `DM` attached to the SNES `must` be a `DMPLEX` at this time.

 89:    Options Database Key:
 90: .  -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`

 92:    Level: advanced

 94: .seealso: `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`, `DMSetNullSpaceConstructor()`
 95: @*/
 96: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
 97: {
 98:   PetscBool monitorIntegral = PETSC_FALSE;

100:   SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx);
101:   if (monitorIntegral) {
102:     Mat          J;
103:     Vec          u;
104:     MatNullSpace nullspace;
105:     const Vec   *nullvecs;
106:     PetscScalar  pintd;

108:     SNESGetSolution(snes, &u);
109:     SNESGetJacobian(snes, &J, NULL, NULL, NULL);
110:     MatGetNullSpace(J, &nullspace);
111:     MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs);
112:     VecDot(nullvecs[0], u, &pintd);
113:     PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd));
114:   }
115:   if (*reason > 0) {
116:     Mat          J;
117:     Vec          u;
118:     MatNullSpace nullspace;
119:     PetscInt     pfield = 1;

121:     SNESGetSolution(snes, &u);
122:     SNESGetJacobian(snes, &J, NULL, NULL, NULL);
123:     MatGetNullSpace(J, &nullspace);
124:     SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx);
125:   }
126:   return 0;
127: }

129: /************************** Interpolation *******************************/

131: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
132: {
133:   PetscBool isPlex;

135:   PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex);
136:   if (isPlex) {
137:     *plex = dm;
138:     PetscObjectReference((PetscObject)dm);
139:   } else {
140:     PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex);
141:     if (!*plex) {
142:       DMConvert(dm, DMPLEX, plex);
143:       PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex);
144:       if (copy) {
145:         DMCopyDMSNES(dm, *plex);
146:         DMCopyAuxiliaryVec(dm, *plex);
147:       }
148:     } else {
149:       PetscObjectReference((PetscObject)*plex);
150:     }
151:   }
152:   return 0;
153: }

155: /*@C
156:   DMInterpolationCreate - Creates a `DMInterpolationInfo` context

158:   Collective

160:   Input Parameter:
161: . comm - the communicator

163:   Output Parameter:
164: . ctx - the context

166:   Level: beginner

168: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
169: @*/
170: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
171: {
173:   PetscNew(ctx);

175:   (*ctx)->comm   = comm;
176:   (*ctx)->dim    = -1;
177:   (*ctx)->nInput = 0;
178:   (*ctx)->points = NULL;
179:   (*ctx)->cells  = NULL;
180:   (*ctx)->n      = -1;
181:   (*ctx)->coords = NULL;
182:   return 0;
183: }

185: /*@C
186:   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context

188:   Not collective

190:   Input Parameters:
191: + ctx - the context
192: - dim - the spatial dimension

194:   Level: intermediate

196: .seealso: `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
197: @*/
198: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
199: {
201:   ctx->dim = dim;
202:   return 0;
203: }

205: /*@C
206:   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context

208:   Not collective

210:   Input Parameter:
211: . ctx - the context

213:   Output Parameter:
214: . dim - the spatial dimension

216:   Level: intermediate

218: .seealso: `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
219: @*/
220: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
221: {
223:   *dim = ctx->dim;
224:   return 0;
225: }

227: /*@C
228:   DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context

230:   Not collective

232:   Input Parameters:
233: + ctx - the context
234: - dof - the number of fields

236:   Level: intermediate

238: .seealso: `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
239: @*/
240: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
241: {
243:   ctx->dof = dof;
244:   return 0;
245: }

247: /*@C
248:   DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context

250:   Not collective

252:   Input Parameter:
253: . ctx - the context

255:   Output Parameter:
256: . dof - the number of fields

258:   Level: intermediate

260: .seealso: DMInterpolationInfo, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
261: @*/
262: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
263: {
265:   *dof = ctx->dof;
266:   return 0;
267: }

269: /*@C
270:   DMInterpolationAddPoints - Add points at which we will interpolate the fields

272:   Not collective

274:   Input Parameters:
275: + ctx    - the context
276: . n      - the number of points
277: - points - the coordinates for each point, an array of size n * dim

279:   Note:
280:   The coordinate information is copied.

282:   Level: intermediate

284: .seealso: `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
285: @*/
286: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
287: {
290:   ctx->nInput = n;

292:   PetscMalloc1(n * ctx->dim, &ctx->points);
293:   PetscArraycpy(ctx->points, points, n * ctx->dim);
294:   return 0;
295: }

297: /*@C
298:   DMInterpolationSetUp - Compute spatial indices for point location during interpolation

300:   Collective

302:   Input Parameters:
303: + ctx - the context
304: . dm  - the `DM` for the function space used for interpolation
305: . redundantPoints - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
306: - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error

308:   Level: intermediate

310: .seealso: DMInterpolationInfo, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
311: @*/
312: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
313: {
314:   MPI_Comm           comm = ctx->comm;
315:   PetscScalar       *a;
316:   PetscInt           p, q, i;
317:   PetscMPIInt        rank, size;
318:   Vec                pointVec;
319:   PetscSF            cellSF;
320:   PetscLayout        layout;
321:   PetscReal         *globalPoints;
322:   PetscScalar       *globalPointsScalar;
323:   const PetscInt    *ranges;
324:   PetscMPIInt       *counts, *displs;
325:   const PetscSFNode *foundCells;
326:   const PetscInt    *foundPoints;
327:   PetscMPIInt       *foundProcs, *globalProcs;
328:   PetscInt           n, N, numFound;

331:   MPI_Comm_size(comm, &size);
332:   MPI_Comm_rank(comm, &rank);
334:   /* Locate points */
335:   n = ctx->nInput;
336:   if (!redundantPoints) {
337:     PetscLayoutCreate(comm, &layout);
338:     PetscLayoutSetBlockSize(layout, 1);
339:     PetscLayoutSetLocalSize(layout, n);
340:     PetscLayoutSetUp(layout);
341:     PetscLayoutGetSize(layout, &N);
342:     /* Communicate all points to all processes */
343:     PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs);
344:     PetscLayoutGetRanges(layout, &ranges);
345:     for (p = 0; p < size; ++p) {
346:       counts[p] = (ranges[p + 1] - ranges[p]) * ctx->dim;
347:       displs[p] = ranges[p] * ctx->dim;
348:     }
349:     MPI_Allgatherv(ctx->points, n * ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
350:   } else {
351:     N            = n;
352:     globalPoints = ctx->points;
353:     counts = displs = NULL;
354:     layout          = NULL;
355:   }
356: #if 0
357:   PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
358:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
359: #else
360:   #if defined(PETSC_USE_COMPLEX)
361:   PetscMalloc1(N * ctx->dim, &globalPointsScalar);
362:   for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
363:   #else
364:   globalPointsScalar = globalPoints;
365:   #endif
366:   VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec);
367:   PetscMalloc2(N, &foundProcs, N, &globalProcs);
368:   for (p = 0; p < N; ++p) foundProcs[p] = size;
369:   cellSF = NULL;
370:   DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF);
371:   PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells);
372: #endif
373:   for (p = 0; p < numFound; ++p) {
374:     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
375:   }
376:   /* Let the lowest rank process own each point */
377:   MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
378:   ctx->n = 0;
379:   for (p = 0; p < N; ++p) {
380:     if (globalProcs[p] == size) {
382:                  (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0));
383:       if (rank == 0) ++ctx->n;
384:     } else if (globalProcs[p] == rank) ++ctx->n;
385:   }
386:   /* Create coordinates vector and array of owned cells */
387:   PetscMalloc1(ctx->n, &ctx->cells);
388:   VecCreate(comm, &ctx->coords);
389:   VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE);
390:   VecSetBlockSize(ctx->coords, ctx->dim);
391:   VecSetType(ctx->coords, VECSTANDARD);
392:   VecGetArray(ctx->coords, &a);
393:   for (p = 0, q = 0, i = 0; p < N; ++p) {
394:     if (globalProcs[p] == rank) {
395:       PetscInt d;

397:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
398:       ctx->cells[q] = foundCells[q].index;
399:       ++q;
400:     }
401:     if (globalProcs[p] == size && rank == 0) {
402:       PetscInt d;

404:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
405:       ctx->cells[q] = -1;
406:       ++q;
407:     }
408:   }
409:   VecRestoreArray(ctx->coords, &a);
410: #if 0
411:   PetscFree3(foundCells,foundProcs,globalProcs);
412: #else
413:   PetscFree2(foundProcs, globalProcs);
414:   PetscSFDestroy(&cellSF);
415:   VecDestroy(&pointVec);
416: #endif
417:   if ((void *)globalPointsScalar != (void *)globalPoints) PetscFree(globalPointsScalar);
418:   if (!redundantPoints) PetscFree3(globalPoints, counts, displs);
419:   PetscLayoutDestroy(&layout);
420:   return 0;
421: }

423: /*@C
424:   DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point

426:   Collective

428:   Input Parameter:
429: . ctx - the context

431:   Output Parameter:
432: . coordinates  - the coordinates of interpolation points

434:   Note:
435:   The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`.
436:   This is a borrowed vector that the user should not destroy.

438:   Level: intermediate

440: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
441: @*/
442: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
443: {
446:   *coordinates = ctx->coords;
447:   return 0;
448: }

450: /*@C
451:   DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values

453:   Collective

455:   Input Parameter:
456: . ctx - the context

458:   Output Parameter:
459: . v  - a vector capable of holding the interpolated field values

461:   Note:
462:   This vector should be returned using `DMInterpolationRestoreVector()`.

464:   Level: intermediate

466: .seealso: `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
467: @*/
468: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
469: {
472:   VecCreate(ctx->comm, v);
473:   VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE);
474:   VecSetBlockSize(*v, ctx->dof);
475:   VecSetType(*v, VECSTANDARD);
476:   return 0;
477: }

479: /*@C
480:   DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values

482:   Collective

484:   Input Parameters:
485: + ctx - the context
486: - v  - a vector capable of holding the interpolated field values

488:   Level: intermediate

490: .seealso: `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
491: @*/
492: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
493: {
496:   VecDestroy(v);
497:   return 0;
498: }

500: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
501: {
502:   PetscReal          v0, J, invJ, detJ;
503:   const PetscInt     dof = ctx->dof;
504:   const PetscScalar *coords;
505:   PetscScalar       *a;
506:   PetscInt           p;

508:   VecGetArrayRead(ctx->coords, &coords);
509:   VecGetArray(v, &a);
510:   for (p = 0; p < ctx->n; ++p) {
511:     PetscInt     c = ctx->cells[p];
512:     PetscScalar *x = NULL;
513:     PetscReal    xir[1];
514:     PetscInt     xSize, comp;

516:     DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ);
518:     xir[0] = invJ * PetscRealPart(coords[p] - v0);
519:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
520:     if (2 * dof == xSize) {
521:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
522:     } else if (dof == xSize) {
523:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
524:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input closure size %" PetscInt_FMT " must be either %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 2 * dof, dof);
525:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
526:   }
527:   VecRestoreArray(v, &a);
528:   VecRestoreArrayRead(ctx->coords, &coords);
529:   return 0;
530: }

532: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
533: {
534:   PetscReal         *v0, *J, *invJ, detJ;
535:   const PetscScalar *coords;
536:   PetscScalar       *a;
537:   PetscInt           p;

539:   PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ);
540:   VecGetArrayRead(ctx->coords, &coords);
541:   VecGetArray(v, &a);
542:   for (p = 0; p < ctx->n; ++p) {
543:     PetscInt     c = ctx->cells[p];
544:     PetscScalar *x = NULL;
545:     PetscReal    xi[4];
546:     PetscInt     d, f, comp;

548:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
550:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
551:     for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];

553:     for (d = 0; d < ctx->dim; ++d) {
554:       xi[d] = 0.0;
555:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
556:       for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[(d + 1) * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
557:     }
558:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
559:   }
560:   VecRestoreArray(v, &a);
561:   VecRestoreArrayRead(ctx->coords, &coords);
562:   PetscFree3(v0, J, invJ);
563:   return 0;
564: }

566: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
567: {
568:   PetscReal         *v0, *J, *invJ, detJ;
569:   const PetscScalar *coords;
570:   PetscScalar       *a;
571:   PetscInt           p;

573:   PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ);
574:   VecGetArrayRead(ctx->coords, &coords);
575:   VecGetArray(v, &a);
576:   for (p = 0; p < ctx->n; ++p) {
577:     PetscInt       c        = ctx->cells[p];
578:     const PetscInt order[3] = {2, 1, 3};
579:     PetscScalar   *x        = NULL;
580:     PetscReal      xi[4];
581:     PetscInt       d, f, comp;

583:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
585:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
586:     for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];

588:     for (d = 0; d < ctx->dim; ++d) {
589:       xi[d] = 0.0;
590:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
591:       for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[order[d] * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
592:     }
593:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
594:   }
595:   VecRestoreArray(v, &a);
596:   VecRestoreArrayRead(ctx->coords, &coords);
597:   PetscFree3(v0, J, invJ);
598:   return 0;
599: }

601: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
602: {
603:   const PetscScalar *vertices = (const PetscScalar *)ctx;
604:   const PetscScalar  x0       = vertices[0];
605:   const PetscScalar  y0       = vertices[1];
606:   const PetscScalar  x1       = vertices[2];
607:   const PetscScalar  y1       = vertices[3];
608:   const PetscScalar  x2       = vertices[4];
609:   const PetscScalar  y2       = vertices[5];
610:   const PetscScalar  x3       = vertices[6];
611:   const PetscScalar  y3       = vertices[7];
612:   const PetscScalar  f_1      = x1 - x0;
613:   const PetscScalar  g_1      = y1 - y0;
614:   const PetscScalar  f_3      = x3 - x0;
615:   const PetscScalar  g_3      = y3 - y0;
616:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
617:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
618:   const PetscScalar *ref;
619:   PetscScalar       *real;

621:   VecGetArrayRead(Xref, &ref);
622:   VecGetArray(Xreal, &real);
623:   {
624:     const PetscScalar p0 = ref[0];
625:     const PetscScalar p1 = ref[1];

627:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
628:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
629:   }
630:   PetscLogFlops(28);
631:   VecRestoreArrayRead(Xref, &ref);
632:   VecRestoreArray(Xreal, &real);
633:   return 0;
634: }

636: #include <petsc/private/dmimpl.h>
637: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
638: {
639:   const PetscScalar *vertices = (const PetscScalar *)ctx;
640:   const PetscScalar  x0       = vertices[0];
641:   const PetscScalar  y0       = vertices[1];
642:   const PetscScalar  x1       = vertices[2];
643:   const PetscScalar  y1       = vertices[3];
644:   const PetscScalar  x2       = vertices[4];
645:   const PetscScalar  y2       = vertices[5];
646:   const PetscScalar  x3       = vertices[6];
647:   const PetscScalar  y3       = vertices[7];
648:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
649:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
650:   const PetscScalar *ref;

652:   VecGetArrayRead(Xref, &ref);
653:   {
654:     const PetscScalar x       = ref[0];
655:     const PetscScalar y       = ref[1];
656:     const PetscInt    rows[2] = {0, 1};
657:     PetscScalar       values[4];

659:     values[0] = (x1 - x0 + f_01 * y) * 0.5;
660:     values[1] = (x3 - x0 + f_01 * x) * 0.5;
661:     values[2] = (y1 - y0 + g_01 * y) * 0.5;
662:     values[3] = (y3 - y0 + g_01 * x) * 0.5;
663:     MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
664:   }
665:   PetscLogFlops(30);
666:   VecRestoreArrayRead(Xref, &ref);
667:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
668:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
669:   return 0;
670: }

672: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
673: {
674:   DM                 dmCoord;
675:   PetscFE            fem = NULL;
676:   SNES               snes;
677:   KSP                ksp;
678:   PC                 pc;
679:   Vec                coordsLocal, r, ref, real;
680:   Mat                J;
681:   PetscTabulation    T = NULL;
682:   const PetscScalar *coords;
683:   PetscScalar       *a;
684:   PetscReal          xir[2] = {0., 0.};
685:   PetscInt           Nf, p;
686:   const PetscInt     dof = ctx->dof;

688:   DMGetNumFields(dm, &Nf);
689:   if (Nf) {
690:     PetscObject  obj;
691:     PetscClassId id;

693:     DMGetField(dm, 0, NULL, &obj);
694:     PetscObjectGetClassId(obj, &id);
695:     if (id == PETSCFE_CLASSID) {
696:       fem = (PetscFE)obj;
697:       PetscFECreateTabulation(fem, 1, 1, xir, 0, &T);
698:     }
699:   }
700:   DMGetCoordinatesLocal(dm, &coordsLocal);
701:   DMGetCoordinateDM(dm, &dmCoord);
702:   SNESCreate(PETSC_COMM_SELF, &snes);
703:   SNESSetOptionsPrefix(snes, "quad_interp_");
704:   VecCreate(PETSC_COMM_SELF, &r);
705:   VecSetSizes(r, 2, 2);
706:   VecSetType(r, dm->vectype);
707:   VecDuplicate(r, &ref);
708:   VecDuplicate(r, &real);
709:   MatCreate(PETSC_COMM_SELF, &J);
710:   MatSetSizes(J, 2, 2, 2, 2);
711:   MatSetType(J, MATSEQDENSE);
712:   MatSetUp(J);
713:   SNESSetFunction(snes, r, QuadMap_Private, NULL);
714:   SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
715:   SNESGetKSP(snes, &ksp);
716:   KSPGetPC(ksp, &pc);
717:   PCSetType(pc, PCLU);
718:   SNESSetFromOptions(snes);

720:   VecGetArrayRead(ctx->coords, &coords);
721:   VecGetArray(v, &a);
722:   for (p = 0; p < ctx->n; ++p) {
723:     PetscScalar *x = NULL, *vertices = NULL;
724:     PetscScalar *xi;
725:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

727:     /* Can make this do all points at once */
728:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
730:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
731:     SNESSetFunction(snes, NULL, NULL, vertices);
732:     SNESSetJacobian(snes, NULL, NULL, NULL, vertices);
733:     VecGetArray(real, &xi);
734:     xi[0] = coords[p * ctx->dim + 0];
735:     xi[1] = coords[p * ctx->dim + 1];
736:     VecRestoreArray(real, &xi);
737:     SNESSolve(snes, real, ref);
738:     VecGetArray(ref, &xi);
739:     xir[0] = PetscRealPart(xi[0]);
740:     xir[1] = PetscRealPart(xi[1]);
741:     if (4 * dof == xSize) {
742:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) * (1 - xir[1]) + x[1 * dof + comp] * xir[0] * (1 - xir[1]) + x[2 * dof + comp] * xir[0] * xir[1] + x[3 * dof + comp] * (1 - xir[0]) * xir[1];
743:     } else if (dof == xSize) {
744:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
745:     } else {
746:       PetscInt d;

749:       xir[0] = 2.0 * xir[0] - 1.0;
750:       xir[1] = 2.0 * xir[1] - 1.0;
751:       PetscFEComputeTabulation(fem, 1, xir, 0, T);
752:       for (comp = 0; comp < dof; ++comp) {
753:         a[p * dof + comp] = 0.0;
754:         for (d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp];
755:       }
756:     }
757:     VecRestoreArray(ref, &xi);
758:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
759:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
760:   }
761:   PetscTabulationDestroy(&T);
762:   VecRestoreArray(v, &a);
763:   VecRestoreArrayRead(ctx->coords, &coords);

765:   SNESDestroy(&snes);
766:   VecDestroy(&r);
767:   VecDestroy(&ref);
768:   VecDestroy(&real);
769:   MatDestroy(&J);
770:   return 0;
771: }

773: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
774: {
775:   const PetscScalar *vertices = (const PetscScalar *)ctx;
776:   const PetscScalar  x0       = vertices[0];
777:   const PetscScalar  y0       = vertices[1];
778:   const PetscScalar  z0       = vertices[2];
779:   const PetscScalar  x1       = vertices[9];
780:   const PetscScalar  y1       = vertices[10];
781:   const PetscScalar  z1       = vertices[11];
782:   const PetscScalar  x2       = vertices[6];
783:   const PetscScalar  y2       = vertices[7];
784:   const PetscScalar  z2       = vertices[8];
785:   const PetscScalar  x3       = vertices[3];
786:   const PetscScalar  y3       = vertices[4];
787:   const PetscScalar  z3       = vertices[5];
788:   const PetscScalar  x4       = vertices[12];
789:   const PetscScalar  y4       = vertices[13];
790:   const PetscScalar  z4       = vertices[14];
791:   const PetscScalar  x5       = vertices[15];
792:   const PetscScalar  y5       = vertices[16];
793:   const PetscScalar  z5       = vertices[17];
794:   const PetscScalar  x6       = vertices[18];
795:   const PetscScalar  y6       = vertices[19];
796:   const PetscScalar  z6       = vertices[20];
797:   const PetscScalar  x7       = vertices[21];
798:   const PetscScalar  y7       = vertices[22];
799:   const PetscScalar  z7       = vertices[23];
800:   const PetscScalar  f_1      = x1 - x0;
801:   const PetscScalar  g_1      = y1 - y0;
802:   const PetscScalar  h_1      = z1 - z0;
803:   const PetscScalar  f_3      = x3 - x0;
804:   const PetscScalar  g_3      = y3 - y0;
805:   const PetscScalar  h_3      = z3 - z0;
806:   const PetscScalar  f_4      = x4 - x0;
807:   const PetscScalar  g_4      = y4 - y0;
808:   const PetscScalar  h_4      = z4 - z0;
809:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
810:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
811:   const PetscScalar  h_01     = z2 - z1 - z3 + z0;
812:   const PetscScalar  f_12     = x7 - x3 - x4 + x0;
813:   const PetscScalar  g_12     = y7 - y3 - y4 + y0;
814:   const PetscScalar  h_12     = z7 - z3 - z4 + z0;
815:   const PetscScalar  f_02     = x5 - x1 - x4 + x0;
816:   const PetscScalar  g_02     = y5 - y1 - y4 + y0;
817:   const PetscScalar  h_02     = z5 - z1 - z4 + z0;
818:   const PetscScalar  f_012    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
819:   const PetscScalar  g_012    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
820:   const PetscScalar  h_012    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
821:   const PetscScalar *ref;
822:   PetscScalar       *real;

824:   VecGetArrayRead(Xref, &ref);
825:   VecGetArray(Xreal, &real);
826:   {
827:     const PetscScalar p0 = ref[0];
828:     const PetscScalar p1 = ref[1];
829:     const PetscScalar p2 = ref[2];

831:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_4 * p2 + f_01 * p0 * p1 + f_12 * p1 * p2 + f_02 * p0 * p2 + f_012 * p0 * p1 * p2;
832:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_4 * p2 + g_01 * p0 * p1 + g_01 * p0 * p1 + g_12 * p1 * p2 + g_02 * p0 * p2 + g_012 * p0 * p1 * p2;
833:     real[2] = z0 + h_1 * p0 + h_3 * p1 + h_4 * p2 + h_01 * p0 * p1 + h_01 * p0 * p1 + h_12 * p1 * p2 + h_02 * p0 * p2 + h_012 * p0 * p1 * p2;
834:   }
835:   PetscLogFlops(114);
836:   VecRestoreArrayRead(Xref, &ref);
837:   VecRestoreArray(Xreal, &real);
838:   return 0;
839: }

841: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
842: {
843:   const PetscScalar *vertices = (const PetscScalar *)ctx;
844:   const PetscScalar  x0       = vertices[0];
845:   const PetscScalar  y0       = vertices[1];
846:   const PetscScalar  z0       = vertices[2];
847:   const PetscScalar  x1       = vertices[9];
848:   const PetscScalar  y1       = vertices[10];
849:   const PetscScalar  z1       = vertices[11];
850:   const PetscScalar  x2       = vertices[6];
851:   const PetscScalar  y2       = vertices[7];
852:   const PetscScalar  z2       = vertices[8];
853:   const PetscScalar  x3       = vertices[3];
854:   const PetscScalar  y3       = vertices[4];
855:   const PetscScalar  z3       = vertices[5];
856:   const PetscScalar  x4       = vertices[12];
857:   const PetscScalar  y4       = vertices[13];
858:   const PetscScalar  z4       = vertices[14];
859:   const PetscScalar  x5       = vertices[15];
860:   const PetscScalar  y5       = vertices[16];
861:   const PetscScalar  z5       = vertices[17];
862:   const PetscScalar  x6       = vertices[18];
863:   const PetscScalar  y6       = vertices[19];
864:   const PetscScalar  z6       = vertices[20];
865:   const PetscScalar  x7       = vertices[21];
866:   const PetscScalar  y7       = vertices[22];
867:   const PetscScalar  z7       = vertices[23];
868:   const PetscScalar  f_xy     = x2 - x1 - x3 + x0;
869:   const PetscScalar  g_xy     = y2 - y1 - y3 + y0;
870:   const PetscScalar  h_xy     = z2 - z1 - z3 + z0;
871:   const PetscScalar  f_yz     = x7 - x3 - x4 + x0;
872:   const PetscScalar  g_yz     = y7 - y3 - y4 + y0;
873:   const PetscScalar  h_yz     = z7 - z3 - z4 + z0;
874:   const PetscScalar  f_xz     = x5 - x1 - x4 + x0;
875:   const PetscScalar  g_xz     = y5 - y1 - y4 + y0;
876:   const PetscScalar  h_xz     = z5 - z1 - z4 + z0;
877:   const PetscScalar  f_xyz    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
878:   const PetscScalar  g_xyz    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
879:   const PetscScalar  h_xyz    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
880:   const PetscScalar *ref;

882:   VecGetArrayRead(Xref, &ref);
883:   {
884:     const PetscScalar x       = ref[0];
885:     const PetscScalar y       = ref[1];
886:     const PetscScalar z       = ref[2];
887:     const PetscInt    rows[3] = {0, 1, 2};
888:     PetscScalar       values[9];

890:     values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
891:     values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
892:     values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
893:     values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
894:     values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
895:     values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
896:     values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
897:     values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
898:     values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;

900:     MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
901:   }
902:   PetscLogFlops(152);
903:   VecRestoreArrayRead(Xref, &ref);
904:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
905:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
906:   return 0;
907: }

909: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
910: {
911:   DM                 dmCoord;
912:   SNES               snes;
913:   KSP                ksp;
914:   PC                 pc;
915:   Vec                coordsLocal, r, ref, real;
916:   Mat                J;
917:   const PetscScalar *coords;
918:   PetscScalar       *a;
919:   PetscInt           p;

921:   DMGetCoordinatesLocal(dm, &coordsLocal);
922:   DMGetCoordinateDM(dm, &dmCoord);
923:   SNESCreate(PETSC_COMM_SELF, &snes);
924:   SNESSetOptionsPrefix(snes, "hex_interp_");
925:   VecCreate(PETSC_COMM_SELF, &r);
926:   VecSetSizes(r, 3, 3);
927:   VecSetType(r, dm->vectype);
928:   VecDuplicate(r, &ref);
929:   VecDuplicate(r, &real);
930:   MatCreate(PETSC_COMM_SELF, &J);
931:   MatSetSizes(J, 3, 3, 3, 3);
932:   MatSetType(J, MATSEQDENSE);
933:   MatSetUp(J);
934:   SNESSetFunction(snes, r, HexMap_Private, NULL);
935:   SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
936:   SNESGetKSP(snes, &ksp);
937:   KSPGetPC(ksp, &pc);
938:   PCSetType(pc, PCLU);
939:   SNESSetFromOptions(snes);

941:   VecGetArrayRead(ctx->coords, &coords);
942:   VecGetArray(v, &a);
943:   for (p = 0; p < ctx->n; ++p) {
944:     PetscScalar *x = NULL, *vertices = NULL;
945:     PetscScalar *xi;
946:     PetscReal    xir[3];
947:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

949:     /* Can make this do all points at once */
950:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
952:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
954:     SNESSetFunction(snes, NULL, NULL, vertices);
955:     SNESSetJacobian(snes, NULL, NULL, NULL, vertices);
956:     VecGetArray(real, &xi);
957:     xi[0] = coords[p * ctx->dim + 0];
958:     xi[1] = coords[p * ctx->dim + 1];
959:     xi[2] = coords[p * ctx->dim + 2];
960:     VecRestoreArray(real, &xi);
961:     SNESSolve(snes, real, ref);
962:     VecGetArray(ref, &xi);
963:     xir[0] = PetscRealPart(xi[0]);
964:     xir[1] = PetscRealPart(xi[1]);
965:     xir[2] = PetscRealPart(xi[2]);
966:     if (8 * ctx->dof == xSize) {
967:       for (comp = 0; comp < ctx->dof; ++comp) {
968:         a[p * ctx->dof + comp] = x[0 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * (1 - xir[2]) + x[3 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * (1 - xir[2]) + x[2 * ctx->dof + comp] * xir[0] * xir[1] * (1 - xir[2]) + x[1 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * (1 - xir[2]) +
969:                                  x[4 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * xir[2] + x[5 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * xir[2] + x[6 * ctx->dof + comp] * xir[0] * xir[1] * xir[2] + x[7 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * xir[2];
970:       }
971:     } else {
972:       for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
973:     }
974:     VecRestoreArray(ref, &xi);
975:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
976:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
977:   }
978:   VecRestoreArray(v, &a);
979:   VecRestoreArrayRead(ctx->coords, &coords);

981:   SNESDestroy(&snes);
982:   VecDestroy(&r);
983:   VecDestroy(&ref);
984:   VecDestroy(&real);
985:   MatDestroy(&J);
986:   return 0;
987: }

989: /*@C
990:   DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.

992:   Input Parameters:
993: + ctx - The `DMInterpolationInfo` context
994: . dm  - The `DM`
995: - x   - The local vector containing the field to be interpolated

997:   Output Parameters:
998: . v   - The vector containing the interpolated values

1000:   Note:
1001:   A suitable v can be obtained using `DMInterpolationGetVector()`.

1003:   Level: beginner

1005: .seealso: `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1006: @*/
1007: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
1008: {
1009:   PetscDS   ds;
1010:   PetscInt  n, p, Nf, field;
1011:   PetscBool useDS = PETSC_FALSE;

1016:   VecGetLocalSize(v, &n);
1018:   if (!n) return 0;
1019:   DMGetDS(dm, &ds);
1020:   if (ds) {
1021:     useDS = PETSC_TRUE;
1022:     PetscDSGetNumFields(ds, &Nf);
1023:     for (field = 0; field < Nf; ++field) {
1024:       PetscObject  obj;
1025:       PetscClassId id;

1027:       PetscDSGetDiscretization(ds, field, &obj);
1028:       PetscObjectGetClassId(obj, &id);
1029:       if (id != PETSCFE_CLASSID) {
1030:         useDS = PETSC_FALSE;
1031:         break;
1032:       }
1033:     }
1034:   }
1035:   if (useDS) {
1036:     const PetscScalar *coords;
1037:     PetscScalar       *interpolant;
1038:     PetscInt           cdim, d;

1040:     DMGetCoordinateDim(dm, &cdim);
1041:     VecGetArrayRead(ctx->coords, &coords);
1042:     VecGetArrayWrite(v, &interpolant);
1043:     for (p = 0; p < ctx->n; ++p) {
1044:       PetscReal    pcoords[3], xi[3];
1045:       PetscScalar *xa   = NULL;
1046:       PetscInt     coff = 0, foff = 0, clSize;

1048:       if (ctx->cells[p] < 0) continue;
1049:       for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
1050:       DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi);
1051:       DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa);
1052:       for (field = 0; field < Nf; ++field) {
1053:         PetscTabulation T;
1054:         PetscFE         fe;

1056:         PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
1057:         PetscFECreateTabulation(fe, 1, 1, xi, 0, &T);
1058:         {
1059:           const PetscReal *basis = T->T[0];
1060:           const PetscInt   Nb    = T->Nb;
1061:           const PetscInt   Nc    = T->Nc;
1062:           PetscInt         f, fc;

1064:           for (fc = 0; fc < Nc; ++fc) {
1065:             interpolant[p * ctx->dof + coff + fc] = 0.0;
1066:             for (f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
1067:           }
1068:           coff += Nc;
1069:           foff += Nb;
1070:         }
1071:         PetscTabulationDestroy(&T);
1072:       }
1073:       DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa);
1076:     }
1077:     VecRestoreArrayRead(ctx->coords, &coords);
1078:     VecRestoreArrayWrite(v, &interpolant);
1079:   } else {
1080:     DMPolytopeType ct;

1082:     /* TODO Check each cell individually */
1083:     DMPlexGetCellType(dm, ctx->cells[0], &ct);
1084:     switch (ct) {
1085:     case DM_POLYTOPE_SEGMENT:
1086:       DMInterpolate_Segment_Private(ctx, dm, x, v);
1087:       break;
1088:     case DM_POLYTOPE_TRIANGLE:
1089:       DMInterpolate_Triangle_Private(ctx, dm, x, v);
1090:       break;
1091:     case DM_POLYTOPE_QUADRILATERAL:
1092:       DMInterpolate_Quad_Private(ctx, dm, x, v);
1093:       break;
1094:     case DM_POLYTOPE_TETRAHEDRON:
1095:       DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
1096:       break;
1097:     case DM_POLYTOPE_HEXAHEDRON:
1098:       DMInterpolate_Hex_Private(ctx, dm, x, v);
1099:       break;
1100:     default:
1101:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
1102:     }
1103:   }
1104:   return 0;
1105: }

1107: /*@C
1108:   DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context

1110:   Collective

1112:   Input Parameter:
1113: . ctx - the context

1115:   Level: beginner

1117: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1118: @*/
1119: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1120: {
1122:   VecDestroy(&(*ctx)->coords);
1123:   PetscFree((*ctx)->points);
1124:   PetscFree((*ctx)->cells);
1125:   PetscFree(*ctx);
1126:   *ctx = NULL;
1127:   return 0;
1128: }

1130: /*@C
1131:   SNESMonitorFields - Monitors the residual for each field separately

1133:   Collective

1135:   Input Parameters:
1136: + snes   - the `SNES` context
1137: . its    - iteration number
1138: . fgnorm - 2-norm of residual
1139: - vf  - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII`

1141:   Note:
1142:   This routine prints the residual norm at each iteration.

1144:   Level: intermediate

1146: .seealso: `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
1147: @*/
1148: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1149: {
1150:   PetscViewer        viewer = vf->viewer;
1151:   Vec                res;
1152:   DM                 dm;
1153:   PetscSection       s;
1154:   const PetscScalar *r;
1155:   PetscReal         *lnorms, *norms;
1156:   PetscInt           numFields, f, pStart, pEnd, p;

1159:   SNESGetFunction(snes, &res, NULL, NULL);
1160:   SNESGetDM(snes, &dm);
1161:   DMGetLocalSection(dm, &s);
1162:   PetscSectionGetNumFields(s, &numFields);
1163:   PetscSectionGetChart(s, &pStart, &pEnd);
1164:   PetscCalloc2(numFields, &lnorms, numFields, &norms);
1165:   VecGetArrayRead(res, &r);
1166:   for (p = pStart; p < pEnd; ++p) {
1167:     for (f = 0; f < numFields; ++f) {
1168:       PetscInt fdof, foff, d;

1170:       PetscSectionGetFieldDof(s, p, f, &fdof);
1171:       PetscSectionGetFieldOffset(s, p, f, &foff);
1172:       for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
1173:     }
1174:   }
1175:   VecRestoreArrayRead(res, &r);
1176:   MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1177:   PetscViewerPushFormat(viewer, vf->format);
1178:   PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel);
1179:   PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm);
1180:   for (f = 0; f < numFields; ++f) {
1181:     if (f > 0) PetscViewerASCIIPrintf(viewer, ", ");
1182:     PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f]));
1183:   }
1184:   PetscViewerASCIIPrintf(viewer, "]\n");
1185:   PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel);
1186:   PetscViewerPopFormat(viewer);
1187:   PetscFree2(lnorms, norms);
1188:   return 0;
1189: }

1191: /********************* Residual Computation **************************/

1193: PetscErrorCode DMPlexGetAllCells_Internal(DM plex, IS *cellIS)
1194: {
1195:   PetscInt depth;

1197:   DMPlexGetDepth(plex, &depth);
1198:   DMGetStratumIS(plex, "dim", depth, cellIS);
1199:   if (!*cellIS) DMGetStratumIS(plex, "depth", depth, cellIS);
1200:   return 0;
1201: }

1203: /*@
1204:   DMPlexSNESComputeResidualFEM - Sums the local residual into vector F from the local input X using pointwise functions specified by the user

1206:   Input Parameters:
1207: + dm - The mesh
1208: . X  - Local solution
1209: - user - The user context

1211:   Output Parameter:
1212: . F  - Local output vector

1214:   Note:
1215:   The residual is summed into F; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in F is intentional.

1217:   Level: developer

1219: .seealso: `DM`, `DMPlexComputeJacobianAction()`
1220: @*/
1221: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1222: {
1223:   DM       plex;
1224:   IS       allcellIS;
1225:   PetscInt Nds, s;

1227:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1228:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1229:   DMGetNumDS(dm, &Nds);
1230:   for (s = 0; s < Nds; ++s) {
1231:     PetscDS      ds;
1232:     IS           cellIS;
1233:     PetscFormKey key;

1235:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1236:     key.value = 0;
1237:     key.field = 0;
1238:     key.part  = 0;
1239:     if (!key.label) {
1240:       PetscObjectReference((PetscObject)allcellIS);
1241:       cellIS = allcellIS;
1242:     } else {
1243:       IS pointIS;

1245:       key.value = 1;
1246:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
1247:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1248:       ISDestroy(&pointIS);
1249:     }
1250:     DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1251:     ISDestroy(&cellIS);
1252:   }
1253:   ISDestroy(&allcellIS);
1254:   DMDestroy(&plex);
1255:   return 0;
1256: }

1258: PetscErrorCode DMSNESComputeResidual(DM dm, Vec X, Vec F, void *user)
1259: {
1260:   DM       plex;
1261:   IS       allcellIS;
1262:   PetscInt Nds, s;

1264:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1265:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1266:   DMGetNumDS(dm, &Nds);
1267:   for (s = 0; s < Nds; ++s) {
1268:     PetscDS ds;
1269:     DMLabel label;
1270:     IS      cellIS;

1272:     DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1273:     {
1274:       PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
1275:       PetscWeakForm     wf;
1276:       PetscInt          Nm = 2, m, Nk = 0, k, kp, off = 0;
1277:       PetscFormKey     *reskeys;

1279:       /* Get unique residual keys */
1280:       for (m = 0; m < Nm; ++m) {
1281:         PetscInt Nkm;
1282:         PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm);
1283:         Nk += Nkm;
1284:       }
1285:       PetscMalloc1(Nk, &reskeys);
1286:       for (m = 0; m < Nm; ++m) PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys);
1288:       PetscFormKeySort(Nk, reskeys);
1289:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1290:         if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
1291:           ++k;
1292:           if (kp != k) reskeys[k] = reskeys[kp];
1293:         }
1294:       }
1295:       Nk = k;

1297:       PetscDSGetWeakForm(ds, &wf);
1298:       for (k = 0; k < Nk; ++k) {
1299:         DMLabel  label = reskeys[k].label;
1300:         PetscInt val   = reskeys[k].value;

1302:         if (!label) {
1303:           PetscObjectReference((PetscObject)allcellIS);
1304:           cellIS = allcellIS;
1305:         } else {
1306:           IS pointIS;

1308:           DMLabelGetStratumIS(label, val, &pointIS);
1309:           ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1310:           ISDestroy(&pointIS);
1311:         }
1312:         DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1313:         ISDestroy(&cellIS);
1314:       }
1315:       PetscFree(reskeys);
1316:     }
1317:   }
1318:   ISDestroy(&allcellIS);
1319:   DMDestroy(&plex);
1320:   return 0;
1321: }

1323: /*@
1324:   DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X

1326:   Input Parameters:
1327: + dm - The mesh
1328: - user - The user context

1330:   Output Parameter:
1331: . X  - Local solution

1333:   Level: developer

1335: .seealso: `DMPLEX`, `DMPlexComputeJacobianAction()`
1336: @*/
1337: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1338: {
1339:   DM plex;

1341:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1342:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
1343:   DMDestroy(&plex);
1344:   return 0;
1345: }

1347: /*@
1348:   DMSNESComputeJacobianAction - Compute the action of the Jacobian J(X) on Y

1350:   Input Parameters:
1351: + dm   - The `DM`
1352: . X    - Local solution vector
1353: . Y    - Local input vector
1354: - user - The user context

1356:   Output Parameter:
1357: . F    - lcoal output vector

1359:   Level: developer

1361:   Notes:
1362:   Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.

1364: .seealso: `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
1365: @*/
1366: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
1367: {
1368:   DM       plex;
1369:   IS       allcellIS;
1370:   PetscInt Nds, s;

1372:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1373:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1374:   DMGetNumDS(dm, &Nds);
1375:   for (s = 0; s < Nds; ++s) {
1376:     PetscDS ds;
1377:     DMLabel label;
1378:     IS      cellIS;

1380:     DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1381:     {
1382:       PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
1383:       PetscWeakForm     wf;
1384:       PetscInt          Nm = 4, m, Nk = 0, k, kp, off = 0;
1385:       PetscFormKey     *jackeys;

1387:       /* Get unique Jacobian keys */
1388:       for (m = 0; m < Nm; ++m) {
1389:         PetscInt Nkm;
1390:         PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm);
1391:         Nk += Nkm;
1392:       }
1393:       PetscMalloc1(Nk, &jackeys);
1394:       for (m = 0; m < Nm; ++m) PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys);
1396:       PetscFormKeySort(Nk, jackeys);
1397:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1398:         if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
1399:           ++k;
1400:           if (kp != k) jackeys[k] = jackeys[kp];
1401:         }
1402:       }
1403:       Nk = k;

1405:       PetscDSGetWeakForm(ds, &wf);
1406:       for (k = 0; k < Nk; ++k) {
1407:         DMLabel  label = jackeys[k].label;
1408:         PetscInt val   = jackeys[k].value;

1410:         if (!label) {
1411:           PetscObjectReference((PetscObject)allcellIS);
1412:           cellIS = allcellIS;
1413:         } else {
1414:           IS pointIS;

1416:           DMLabelGetStratumIS(label, val, &pointIS);
1417:           ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1418:           ISDestroy(&pointIS);
1419:         }
1420:         DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user);
1421:         ISDestroy(&cellIS);
1422:       }
1423:       PetscFree(jackeys);
1424:     }
1425:   }
1426:   ISDestroy(&allcellIS);
1427:   DMDestroy(&plex);
1428:   return 0;
1429: }

1431: /*@
1432:   DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.

1434:   Input Parameters:
1435: + dm - The mesh
1436: . X  - Local input vector
1437: - user - The user context

1439:   Output Parameter:
1440: . Jac  - Jacobian matrix

1442:   Note:
1443:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1444:   like a GPU, or vectorize on a multicore machine.

1446:   Level: developer

1448: .seealso: `DMPLEX`, `Mat`
1449: @*/
1450: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
1451: {
1452:   DM        plex;
1453:   IS        allcellIS;
1454:   PetscBool hasJac, hasPrec;
1455:   PetscInt  Nds, s;

1457:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1458:   DMPlexGetAllCells_Internal(plex, &allcellIS);
1459:   DMGetNumDS(dm, &Nds);
1460:   for (s = 0; s < Nds; ++s) {
1461:     PetscDS      ds;
1462:     IS           cellIS;
1463:     PetscFormKey key;

1465:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1466:     key.value = 0;
1467:     key.field = 0;
1468:     key.part  = 0;
1469:     if (!key.label) {
1470:       PetscObjectReference((PetscObject)allcellIS);
1471:       cellIS = allcellIS;
1472:     } else {
1473:       IS pointIS;

1475:       key.value = 1;
1476:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
1477:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1478:       ISDestroy(&pointIS);
1479:     }
1480:     if (!s) {
1481:       PetscDSHasJacobian(ds, &hasJac);
1482:       PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1483:       if (hasJac && hasPrec) MatZeroEntries(Jac);
1484:       MatZeroEntries(JacP);
1485:     }
1486:     DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user);
1487:     ISDestroy(&cellIS);
1488:   }
1489:   ISDestroy(&allcellIS);
1490:   DMDestroy(&plex);
1491:   return 0;
1492: }

1494: struct _DMSNESJacobianMFCtx {
1495:   DM    dm;
1496:   Vec   X;
1497:   void *ctx;
1498: };

1500: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
1501: {
1502:   struct _DMSNESJacobianMFCtx *ctx;

1504:   MatShellGetContext(A, &ctx);
1505:   MatShellSetContext(A, NULL);
1506:   DMDestroy(&ctx->dm);
1507:   VecDestroy(&ctx->X);
1508:   PetscFree(ctx);
1509:   return 0;
1510: }

1512: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
1513: {
1514:   struct _DMSNESJacobianMFCtx *ctx;

1516:   MatShellGetContext(A, &ctx);
1517:   DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx);
1518:   return 0;
1519: }

1521: /*@
1522:   DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free

1524:   Collective

1526:   Input Parameters:
1527: + dm   - The `DM`
1528: . X    - The evaluation point for the Jacobian
1529: - user - A user context, or NULL

1531:   Output Parameter:
1532: . J    - The `Mat`

1534:   Level: advanced

1536:   Note:
1537:   Vec X is kept in `Mat` J, so updating X then updates the evaluation point.

1539: .seealso: `DM`, `DMSNESComputeJacobianAction()`
1540: @*/
1541: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
1542: {
1543:   struct _DMSNESJacobianMFCtx *ctx;
1544:   PetscInt                     n, N;

1546:   MatCreate(PetscObjectComm((PetscObject)dm), J);
1547:   MatSetType(*J, MATSHELL);
1548:   VecGetLocalSize(X, &n);
1549:   VecGetSize(X, &N);
1550:   MatSetSizes(*J, n, n, N, N);
1551:   PetscObjectReference((PetscObject)dm);
1552:   PetscObjectReference((PetscObject)X);
1553:   PetscMalloc1(1, &ctx);
1554:   ctx->dm  = dm;
1555:   ctx->X   = X;
1556:   ctx->ctx = user;
1557:   MatShellSetContext(*J, ctx);
1558:   MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private);
1559:   MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private);
1560:   return 0;
1561: }

1563: /*
1564:      MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context.

1566:    Input Parameters:
1567: +     X - SNES linearization point
1568: .     ovl - index set of overlapping subdomains

1570:    Output Parameter:
1571: .     J - unassembled (Neumann) local matrix

1573:    Level: intermediate

1575: .seealso: `DMCreateNeumannOverlap()`, `MATIS`, `PCHPDDMSetAuxiliaryMat()`
1576: */
1577: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1578: {
1579:   SNES   snes;
1580:   Mat    pJ;
1581:   DM     ovldm, origdm;
1582:   DMSNES sdm;
1583:   PetscErrorCode (*bfun)(DM, Vec, void *);
1584:   PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
1585:   void *bctx, *jctx;

1587:   PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ);
1589:   PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm);
1591:   MatGetDM(pJ, &ovldm);
1592:   DMSNESGetBoundaryLocal(origdm, &bfun, &bctx);
1593:   DMSNESSetBoundaryLocal(ovldm, bfun, bctx);
1594:   DMSNESGetJacobianLocal(origdm, &jfun, &jctx);
1595:   DMSNESSetJacobianLocal(ovldm, jfun, jctx);
1596:   PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes);
1597:   if (!snes) {
1598:     SNESCreate(PetscObjectComm((PetscObject)ovl), &snes);
1599:     SNESSetDM(snes, ovldm);
1600:     PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes);
1601:     PetscObjectDereference((PetscObject)snes);
1602:   }
1603:   DMGetDMSNES(ovldm, &sdm);
1604:   VecLockReadPush(X);
1605:   {
1606:     void *ctx;
1607:     PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
1608:     DMSNESGetJacobian(ovldm, &J, &ctx);
1609:     PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
1610:   }
1611:   VecLockReadPop(X);
1612:   /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1613:   {
1614:     Mat locpJ;

1616:     MatISGetLocalMat(pJ, &locpJ);
1617:     MatCopy(locpJ, J, SAME_NONZERO_PATTERN);
1618:   }
1619:   return 0;
1620: }

1622: /*@
1623:   DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, residual, and Jacobian.

1625:   Input Parameters:
1626: + dm - The `DM` object
1627: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see `PetscDSAddBoundary()`)
1628: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see `PetscDSSetResidual()`)
1629: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see `PetscDSSetJacobian()`)

1631:   Level: developer

1633: .seealso: `DMPLEX`, `SNES`
1634: @*/
1635: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1636: {
1637:   DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, boundaryctx);
1638:   DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, residualctx);
1639:   DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, jacobianctx);
1640:   PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex);
1641:   return 0;
1642: }

1644: /*@C
1645:   DMSNESCheckDiscretization - Check the discretization error of the exact solution

1647:   Input Parameters:
1648: + snes - the `SNES` object
1649: . dm   - the `DM`
1650: . t    - the time
1651: . u    - a `DM` vector
1652: - tol  - A tolerance for the check, or -1 to print the results instead

1654:   Output Parameters:
1655: . error - An array which holds the discretization error in each field, or NULL

1657:   Note:
1658:   The user must call `PetscDSSetExactSolution()` beforehand

1660:   Level: developer

1662: .seealso: `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`, `PetscDSSetExactSolution()`
1663: @*/
1664: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1665: {
1666:   PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1667:   void     **ectxs;
1668:   PetscReal *err;
1669:   MPI_Comm   comm;
1670:   PetscInt   Nf, f;


1677:   DMComputeExactSolution(dm, t, u, NULL);
1678:   VecViewFromOptions(u, NULL, "-vec_view");

1680:   PetscObjectGetComm((PetscObject)snes, &comm);
1681:   DMGetNumFields(dm, &Nf);
1682:   PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);
1683:   {
1684:     PetscInt Nds, s;

1686:     DMGetNumDS(dm, &Nds);
1687:     for (s = 0; s < Nds; ++s) {
1688:       PetscDS         ds;
1689:       DMLabel         label;
1690:       IS              fieldIS;
1691:       const PetscInt *fields;
1692:       PetscInt        dsNf, f;

1694:       DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);
1695:       PetscDSGetNumFields(ds, &dsNf);
1696:       ISGetIndices(fieldIS, &fields);
1697:       for (f = 0; f < dsNf; ++f) {
1698:         const PetscInt field = fields[f];
1699:         PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);
1700:       }
1701:       ISRestoreIndices(fieldIS, &fields);
1702:     }
1703:   }
1704:   if (Nf > 1) {
1705:     DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err);
1706:     if (tol >= 0.0) {
1708:     } else if (error) {
1709:       for (f = 0; f < Nf; ++f) error[f] = err[f];
1710:     } else {
1711:       PetscPrintf(comm, "L_2 Error: [");
1712:       for (f = 0; f < Nf; ++f) {
1713:         if (f) PetscPrintf(comm, ", ");
1714:         PetscPrintf(comm, "%g", (double)err[f]);
1715:       }
1716:       PetscPrintf(comm, "]\n");
1717:     }
1718:   } else {
1719:     DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]);
1720:     if (tol >= 0.0) {
1722:     } else if (error) {
1723:       error[0] = err[0];
1724:     } else {
1725:       PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]);
1726:     }
1727:   }
1728:   PetscFree3(exacts, ectxs, err);
1729:   return 0;
1730: }

1732: /*@C
1733:   DMSNESCheckResidual - Check the residual of the exact solution

1735:   Input Parameters:
1736: + snes - the `SNES` object
1737: . dm   - the `DM`
1738: . u    - a `DM` vector
1739: - tol  - A tolerance for the check, or -1 to print the results instead

1741:   Output Parameter:
1742: . residual - The residual norm of the exact solution, or NULL

1744:   Level: developer

1746: .seealso: `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
1747: @*/
1748: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1749: {
1750:   MPI_Comm  comm;
1751:   Vec       r;
1752:   PetscReal res;

1758:   PetscObjectGetComm((PetscObject)snes, &comm);
1759:   DMComputeExactSolution(dm, 0.0, u, NULL);
1760:   VecDuplicate(u, &r);
1761:   SNESComputeFunction(snes, u, r);
1762:   VecNorm(r, NORM_2, &res);
1763:   if (tol >= 0.0) {
1765:   } else if (residual) {
1766:     *residual = res;
1767:   } else {
1768:     PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
1769:     VecChop(r, 1.0e-10);
1770:     PetscObjectSetName((PetscObject)r, "Initial Residual");
1771:     PetscObjectSetOptionsPrefix((PetscObject)r, "res_");
1772:     VecViewFromOptions(r, NULL, "-vec_view");
1773:   }
1774:   VecDestroy(&r);
1775:   return 0;
1776: }

1778: /*@C
1779:   DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test

1781:   Input Parameters:
1782: + snes - the `SNES` object
1783: . dm   - the `DM`
1784: . u    - a `DM` vector
1785: - tol  - A tolerance for the check, or -1 to print the results instead

1787:   Output Parameters:
1788: + isLinear - Flag indicaing that the function looks linear, or NULL
1789: - convRate - The rate of convergence of the linear model, or NULL

1791:   Level: developer

1793: .seealso: `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
1794: @*/
1795: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1796: {
1797:   MPI_Comm     comm;
1798:   PetscDS      ds;
1799:   Mat          J, M;
1800:   MatNullSpace nullspace;
1801:   PetscReal    slope, intercept;
1802:   PetscBool    hasJac, hasPrec, isLin = PETSC_FALSE;

1809:   PetscObjectGetComm((PetscObject)snes, &comm);
1810:   DMComputeExactSolution(dm, 0.0, u, NULL);
1811:   /* Create and view matrices */
1812:   DMCreateMatrix(dm, &J);
1813:   DMGetDS(dm, &ds);
1814:   PetscDSHasJacobian(ds, &hasJac);
1815:   PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1816:   if (hasJac && hasPrec) {
1817:     DMCreateMatrix(dm, &M);
1818:     SNESComputeJacobian(snes, u, J, M);
1819:     PetscObjectSetName((PetscObject)M, "Preconditioning Matrix");
1820:     PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_");
1821:     MatViewFromOptions(M, NULL, "-mat_view");
1822:     MatDestroy(&M);
1823:   } else {
1824:     SNESComputeJacobian(snes, u, J, J);
1825:   }
1826:   PetscObjectSetName((PetscObject)J, "Jacobian");
1827:   PetscObjectSetOptionsPrefix((PetscObject)J, "jac_");
1828:   MatViewFromOptions(J, NULL, "-mat_view");
1829:   /* Check nullspace */
1830:   MatGetNullSpace(J, &nullspace);
1831:   if (nullspace) {
1832:     PetscBool isNull;
1833:     MatNullSpaceTest(nullspace, J, &isNull);
1835:   }
1836:   /* Taylor test */
1837:   {
1838:     PetscRandom rand;
1839:     Vec         du, uhat, r, rhat, df;
1840:     PetscReal   h;
1841:     PetscReal  *es, *hs, *errors;
1842:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1843:     PetscInt    Nv, v;

1845:     /* Choose a perturbation direction */
1846:     PetscRandomCreate(comm, &rand);
1847:     VecDuplicate(u, &du);
1848:     VecSetRandom(du, rand);
1849:     PetscRandomDestroy(&rand);
1850:     VecDuplicate(u, &df);
1851:     MatMult(J, du, df);
1852:     /* Evaluate residual at u, F(u), save in vector r */
1853:     VecDuplicate(u, &r);
1854:     SNESComputeFunction(snes, u, r);
1855:     /* Look at the convergence of our Taylor approximation as we approach u */
1856:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
1857:       ;
1858:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
1859:     VecDuplicate(u, &uhat);
1860:     VecDuplicate(u, &rhat);
1861:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1862:       VecWAXPY(uhat, h, du, u);
1863:       /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1864:       SNESComputeFunction(snes, uhat, rhat);
1865:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
1866:       VecNorm(rhat, NORM_2, &errors[Nv]);

1868:       es[Nv] = PetscLog10Real(errors[Nv]);
1869:       hs[Nv] = PetscLog10Real(h);
1870:     }
1871:     VecDestroy(&uhat);
1872:     VecDestroy(&rhat);
1873:     VecDestroy(&df);
1874:     VecDestroy(&r);
1875:     VecDestroy(&du);
1876:     for (v = 0; v < Nv; ++v) {
1877:       if ((tol >= 0) && (errors[v] > tol)) break;
1878:       else if (errors[v] > PETSC_SMALL) break;
1879:     }
1880:     if (v == Nv) isLin = PETSC_TRUE;
1881:     PetscLinearRegression(Nv, hs, es, &slope, &intercept);
1882:     PetscFree3(es, hs, errors);
1883:     /* Slope should be about 2 */
1884:     if (tol >= 0) {
1886:     } else if (isLinear || convRate) {
1887:       if (isLinear) *isLinear = isLin;
1888:       if (convRate) *convRate = slope;
1889:     } else {
1890:       if (!isLin) PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope);
1891:       else PetscPrintf(comm, "Function appears to be linear\n");
1892:     }
1893:   }
1894:   MatDestroy(&J);
1895:   return 0;
1896: }

1898: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1899: {
1900:   DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL);
1901:   DMSNESCheckResidual(snes, dm, u, -1.0, NULL);
1902:   DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL);
1903:   return 0;
1904: }

1906: /*@C
1907:   DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information

1909:   Input Parameters:
1910: + snes - the `SNES` object
1911: - u    - representative `SNES` vector

1913:   Note:
1914:   The user must call `PetscDSSetExactSolution()` beforehand

1916:   Level: developer
1917: @*/
1918: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1919: {
1920:   DM        dm;
1921:   Vec       sol;
1922:   PetscBool check;

1924:   PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check);
1925:   if (!check) return 0;
1926:   SNESGetDM(snes, &dm);
1927:   VecDuplicate(u, &sol);
1928:   SNESSetSolution(snes, sol);
1929:   DMSNESCheck_Internal(snes, dm, sol);
1930:   VecDestroy(&sol);
1931:   return 0;
1932: }