Actual source code: dmplexsnes.c

petsc-master 2020-10-20
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/snesimpl.h>
  3: #include <petscds.h>
  4: #include <petscblaslapack.h>
  5: #include <petsc/private/petscimpl.h>
  6: #include <petsc/private/petscfeimpl.h>

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

 16: /*
 17:   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.

 19:   Collective on SNES

 21:   Input Parameters:
 22: + snes      - The SNES
 23: . pfield    - The field number for pressure
 24: . nullspace - The pressure nullspace
 25: . u         - The solution vector
 26: - ctx       - An optional user context

 28:   Output Parameter:
 29: . u         - The solution with a continuum pressure integral of zero

 31:   Notes:
 32:   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.

 34:   Level: developer

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

 49:   PetscObjectGetComm((PetscObject) snes, &comm);
 50:   SNESGetDM(snes, &dm);
 51:   if (!dm) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM");
 52:   if (!nullspace) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace");
 53:   DMGetDS(dm, &ds);
 54:   PetscDSSetObjective(ds, pfield, pressure_Private);
 55:   MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs);
 56:   if (Nv != 1) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %D", Nv);
 57:   VecDot(nullvecs[0], u, &pintd);
 58:   if (PetscAbsScalar(pintd) > PETSC_SMALL) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g\n", (double) PetscRealPart(pintd));
 59:   PetscDSGetNumFields(ds, &Nf);
 60:   PetscMalloc2(Nf, &intc, Nf, &intn);
 61:   DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx);
 62:   DMPlexComputeIntegralFEM(dm, u, intc, ctx);
 63:   VecAXPY(u, -intc[pfield]/intn[pfield], nullvecs[0]);
 64: #if defined (PETSC_USE_DEBUG)
 65:   DMPlexComputeIntegralFEM(dm, u, intc, ctx);
 66:   if (PetscAbsScalar(intc[pfield]) > PETSC_SMALL) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g\n", (double) PetscRealPart(intc[pfield]));
 67: #endif
 68:   PetscFree2(intc, intn);
 69:   return(0);
 70: }

 72: /*@C
 73:    SNESConvergedCorrectPressure - Convergence test that adds 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. The convergence test itself just mimics SNESConvergedDefault().

 75:    Logically Collective on SNES

 77:    Input Parameters:
 78: +  snes - the SNES context
 79: .  it - the iteration (0 indicates before any Newton steps)
 80: .  xnorm - 2-norm of current iterate
 81: .  snorm - 2-norm of current step
 82: .  fnorm - 2-norm of function at current iterate
 83: -  ctx   - Optional user context

 85:    Output Parameter:
 86: .  reason  - SNES_CONVERGED_ITERATING, SNES_CONVERGED_ITS, or SNES_DIVERGED_FNORM_NAN

 88:    Notes:
 89:    In order to use this monitor, you must setup several PETSc structures. First fields must be added to the DM, and a PetscDS must be created with discretizations of those fields. We currently assume that the pressure field has index 1. The pressure field must have a nullspace, likely created using the DMSetNullSpaceConstructor() interface. Last we must be able to integrate the pressure over the domain, so the DM attached to the SNES must be a Plex at this time.

 91:    Level: advanced

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

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

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

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

130: /************************** Interpolation *******************************/

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

138:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
139:   if (isPlex) {
140:     *plex = dm;
141:     PetscObjectReference((PetscObject) dm);
142:   } else {
143:     PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
144:     if (!*plex) {
145:       DMConvert(dm,DMPLEX,plex);
146:       PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
147:       if (copy) {
148:         PetscInt    i;
149:         PetscObject obj;
150:         const char *comps[3] = {"A","dmAux","dmCh"};

152:         DMCopyDMSNES(dm, *plex);
153:         for (i = 0; i < 3; i++) {
154:           PetscObjectQuery((PetscObject) dm, comps[i], &obj);
155:           PetscObjectCompose((PetscObject) *plex, comps[i], obj);
156:         }
157:       }
158:     } else {
159:       PetscObjectReference((PetscObject) *plex);
160:     }
161:   }
162:   return(0);
163: }

165: /*@C
166:   DMInterpolationCreate - Creates a DMInterpolationInfo context

168:   Collective

170:   Input Parameter:
171: . comm - the communicator

173:   Output Parameter:
174: . ctx - the context

176:   Level: beginner

178: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationDestroy()
179: @*/
180: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
181: {

186:   PetscNew(ctx);

188:   (*ctx)->comm   = comm;
189:   (*ctx)->dim    = -1;
190:   (*ctx)->nInput = 0;
191:   (*ctx)->points = NULL;
192:   (*ctx)->cells  = NULL;
193:   (*ctx)->n      = -1;
194:   (*ctx)->coords = NULL;
195:   return(0);
196: }

198: /*@C
199:   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context

201:   Not collective

203:   Input Parameters:
204: + ctx - the context
205: - dim - the spatial dimension

207:   Level: intermediate

209: .seealso: DMInterpolationGetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
210: @*/
211: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
212: {
214:   if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %D", dim);
215:   ctx->dim = dim;
216:   return(0);
217: }

219: /*@C
220:   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context

222:   Not collective

224:   Input Parameter:
225: . ctx - the context

227:   Output Parameter:
228: . dim - the spatial dimension

230:   Level: intermediate

232: .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
233: @*/
234: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
235: {
238:   *dim = ctx->dim;
239:   return(0);
240: }

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

245:   Not collective

247:   Input Parameters:
248: + ctx - the context
249: - dof - the number of fields

251:   Level: intermediate

253: .seealso: DMInterpolationGetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
254: @*/
255: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
256: {
258:   if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %D", dof);
259:   ctx->dof = dof;
260:   return(0);
261: }

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

266:   Not collective

268:   Input Parameter:
269: . ctx - the context

271:   Output Parameter:
272: . dof - the number of fields

274:   Level: intermediate

276: .seealso: DMInterpolationSetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
277: @*/
278: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
279: {
282:   *dof = ctx->dof;
283:   return(0);
284: }

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

289:   Not collective

291:   Input Parameters:
292: + ctx    - the context
293: . n      - the number of points
294: - points - the coordinates for each point, an array of size n * dim

296:   Note: The coordinate information is copied.

298:   Level: intermediate

300: .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationCreate()
301: @*/
302: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
303: {

307:   if (ctx->dim < 0) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
308:   if (ctx->points)  SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
309:   ctx->nInput = n;

311:   PetscMalloc1(n*ctx->dim, &ctx->points);
312:   PetscArraycpy(ctx->points, points, n*ctx->dim);
313:   return(0);
314: }

316: /*@C
317:   DMInterpolationSetUp - Computea spatial indices that add in point location during interpolation

319:   Collective on ctx

321:   Input Parameters:
322: + ctx - the context
323: . dm  - the DM for the function space used for interpolation
324: - redundantPoints - If PETSC_TRUE, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.

326:   Level: intermediate

328: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
329: @*/
330: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
331: {
332:   MPI_Comm          comm = ctx->comm;
333:   PetscScalar       *a;
334:   PetscInt          p, q, i;
335:   PetscMPIInt       rank, size;
336:   PetscErrorCode    ierr;
337:   Vec               pointVec;
338:   PetscSF           cellSF;
339:   PetscLayout       layout;
340:   PetscReal         *globalPoints;
341:   PetscScalar       *globalPointsScalar;
342:   const PetscInt    *ranges;
343:   PetscMPIInt       *counts, *displs;
344:   const PetscSFNode *foundCells;
345:   const PetscInt    *foundPoints;
346:   PetscMPIInt       *foundProcs, *globalProcs;
347:   PetscInt          n, N, numFound;

351:   MPI_Comm_size(comm, &size);
352:   MPI_Comm_rank(comm, &rank);
353:   if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
354:   /* Locate points */
355:   n = ctx->nInput;
356:   if (!redundantPoints) {
357:     PetscLayoutCreate(comm, &layout);
358:     PetscLayoutSetBlockSize(layout, 1);
359:     PetscLayoutSetLocalSize(layout, n);
360:     PetscLayoutSetUp(layout);
361:     PetscLayoutGetSize(layout, &N);
362:     /* Communicate all points to all processes */
363:     PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
364:     PetscLayoutGetRanges(layout, &ranges);
365:     for (p = 0; p < size; ++p) {
366:       counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
367:       displs[p] = ranges[p]*ctx->dim;
368:     }
369:     MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
370:   } else {
371:     N = n;
372:     globalPoints = ctx->points;
373:     counts = displs = NULL;
374:     layout = NULL;
375:   }
376: #if 0
377:   PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
378:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
379: #else
380: #if defined(PETSC_USE_COMPLEX)
381:   PetscMalloc1(N*ctx->dim,&globalPointsScalar);
382:   for (i=0; i<N*ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
383: #else
384:   globalPointsScalar = globalPoints;
385: #endif
386:   VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
387:   PetscMalloc2(N,&foundProcs,N,&globalProcs);
388:   for (p = 0; p < N; ++p) {foundProcs[p] = size;}
389:   cellSF = NULL;
390:   DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF);
391:   PetscSFGetGraph(cellSF,NULL,&numFound,&foundPoints,&foundCells);
392: #endif
393:   for (p = 0; p < numFound; ++p) {
394:     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
395:   }
396:   /* Let the lowest rank process own each point */
397:   MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
398:   ctx->n = 0;
399:   for (p = 0; p < N; ++p) {
400:     if (globalProcs[p] == size) SETERRQ4(comm, PETSC_ERR_PLIB, "Point %d: %g %g %g not located in mesh", p, (double)globalPoints[p*ctx->dim+0], (double)(ctx->dim > 1 ? globalPoints[p*ctx->dim+1] : 0.0), (double)(ctx->dim > 2 ? globalPoints[p*ctx->dim+2] : 0.0));
401:     else if (globalProcs[p] == rank) ctx->n++;
402:   }
403:   /* Create coordinates vector and array of owned cells */
404:   PetscMalloc1(ctx->n, &ctx->cells);
405:   VecCreate(comm, &ctx->coords);
406:   VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
407:   VecSetBlockSize(ctx->coords, ctx->dim);
408:   VecSetType(ctx->coords,VECSTANDARD);
409:   VecGetArray(ctx->coords, &a);
410:   for (p = 0, q = 0, i = 0; p < N; ++p) {
411:     if (globalProcs[p] == rank) {
412:       PetscInt d;

414:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
415:       ctx->cells[q] = foundCells[q].index;
416:       ++q;
417:     }
418:   }
419:   VecRestoreArray(ctx->coords, &a);
420: #if 0
421:   PetscFree3(foundCells,foundProcs,globalProcs);
422: #else
423:   PetscFree2(foundProcs,globalProcs);
424:   PetscSFDestroy(&cellSF);
425:   VecDestroy(&pointVec);
426: #endif
427:   if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
428:   if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
429:   PetscLayoutDestroy(&layout);
430:   return(0);
431: }

433: /*@C
434:   DMInterpolationGetCoordinates - Gets a Vec with the coordinates of each interpolation point

436:   Collective on ctx

438:   Input Parameter:
439: . ctx - the context

441:   Output Parameter:
442: . coordinates  - the coordinates of interpolation points

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

446:   Level: intermediate

448: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
449: @*/
450: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
451: {
454:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
455:   *coordinates = ctx->coords;
456:   return(0);
457: }

459: /*@C
460:   DMInterpolationGetVector - Gets a Vec which can hold all the interpolated field values

462:   Collective on ctx

464:   Input Parameter:
465: . ctx - the context

467:   Output Parameter:
468: . v  - a vector capable of holding the interpolated field values

470:   Note: This vector should be returned using DMInterpolationRestoreVector().

472:   Level: intermediate

474: .seealso: DMInterpolationRestoreVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
475: @*/
476: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
477: {

482:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
483:   VecCreate(ctx->comm, v);
484:   VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
485:   VecSetBlockSize(*v, ctx->dof);
486:   VecSetType(*v,VECSTANDARD);
487:   return(0);
488: }

490: /*@C
491:   DMInterpolationRestoreVector - Returns a Vec which can hold all the interpolated field values

493:   Collective on ctx

495:   Input Parameters:
496: + ctx - the context
497: - v  - a vector capable of holding the interpolated field values

499:   Level: intermediate

501: .seealso: DMInterpolationGetVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
502: @*/
503: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
504: {

509:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
510:   VecDestroy(v);
511:   return(0);
512: }

514: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
515: {
516:   PetscReal      *v0, *J, *invJ, detJ;
517:   const PetscScalar *coords;
518:   PetscScalar    *a;
519:   PetscInt       p;

523:   PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
524:   VecGetArrayRead(ctx->coords, &coords);
525:   VecGetArray(v, &a);
526:   for (p = 0; p < ctx->n; ++p) {
527:     PetscInt     c = ctx->cells[p];
528:     PetscScalar *x = NULL;
529:     PetscReal    xi[4];
530:     PetscInt     d, f, comp;

532:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
533:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double)detJ, c);
534:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
535:     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];

537:     for (d = 0; d < ctx->dim; ++d) {
538:       xi[d] = 0.0;
539:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
540:       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];
541:     }
542:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
543:   }
544:   VecRestoreArray(v, &a);
545:   VecRestoreArrayRead(ctx->coords, &coords);
546:   PetscFree3(v0, J, invJ);
547:   return(0);
548: }

550: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
551: {
552:   PetscReal      *v0, *J, *invJ, detJ;
553:   const PetscScalar *coords;
554:   PetscScalar    *a;
555:   PetscInt       p;

559:   PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
560:   VecGetArrayRead(ctx->coords, &coords);
561:   VecGetArray(v, &a);
562:   for (p = 0; p < ctx->n; ++p) {
563:     PetscInt       c = ctx->cells[p];
564:     const PetscInt order[3] = {2, 1, 3};
565:     PetscScalar   *x = NULL;
566:     PetscReal      xi[4];
567:     PetscInt       d, f, comp;

569:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
570:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double)detJ, c);
571:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
572:     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];

574:     for (d = 0; d < ctx->dim; ++d) {
575:       xi[d] = 0.0;
576:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
577:       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];
578:     }
579:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
580:   }
581:   VecRestoreArray(v, &a);
582:   VecRestoreArrayRead(ctx->coords, &coords);
583:   PetscFree3(v0, J, invJ);
584:   return(0);
585: }

587: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
588: {
589:   const PetscScalar *vertices = (const PetscScalar*) ctx;
590:   const PetscScalar x0        = vertices[0];
591:   const PetscScalar y0        = vertices[1];
592:   const PetscScalar x1        = vertices[2];
593:   const PetscScalar y1        = vertices[3];
594:   const PetscScalar x2        = vertices[4];
595:   const PetscScalar y2        = vertices[5];
596:   const PetscScalar x3        = vertices[6];
597:   const PetscScalar y3        = vertices[7];
598:   const PetscScalar f_1       = x1 - x0;
599:   const PetscScalar g_1       = y1 - y0;
600:   const PetscScalar f_3       = x3 - x0;
601:   const PetscScalar g_3       = y3 - y0;
602:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
603:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
604:   const PetscScalar *ref;
605:   PetscScalar       *real;
606:   PetscErrorCode    ierr;

609:   VecGetArrayRead(Xref,  &ref);
610:   VecGetArray(Xreal, &real);
611:   {
612:     const PetscScalar p0 = ref[0];
613:     const PetscScalar p1 = ref[1];

615:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
616:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
617:   }
618:   PetscLogFlops(28);
619:   VecRestoreArrayRead(Xref,  &ref);
620:   VecRestoreArray(Xreal, &real);
621:   return(0);
622: }

624: #include <petsc/private/dmimpl.h>
625: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
626: {
627:   const PetscScalar *vertices = (const PetscScalar*) ctx;
628:   const PetscScalar x0        = vertices[0];
629:   const PetscScalar y0        = vertices[1];
630:   const PetscScalar x1        = vertices[2];
631:   const PetscScalar y1        = vertices[3];
632:   const PetscScalar x2        = vertices[4];
633:   const PetscScalar y2        = vertices[5];
634:   const PetscScalar x3        = vertices[6];
635:   const PetscScalar y3        = vertices[7];
636:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
637:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
638:   const PetscScalar *ref;
639:   PetscErrorCode    ierr;

642:   VecGetArrayRead(Xref,  &ref);
643:   {
644:     const PetscScalar x       = ref[0];
645:     const PetscScalar y       = ref[1];
646:     const PetscInt    rows[2] = {0, 1};
647:     PetscScalar       values[4];

649:     values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
650:     values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
651:     MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
652:   }
653:   PetscLogFlops(30);
654:   VecRestoreArrayRead(Xref,  &ref);
655:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
656:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
657:   return(0);
658: }

660: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
661: {
662:   DM             dmCoord;
663:   PetscFE        fem = NULL;
664:   SNES           snes;
665:   KSP            ksp;
666:   PC             pc;
667:   Vec            coordsLocal, r, ref, real;
668:   Mat            J;
669:   PetscTabulation    T;
670:   const PetscScalar *coords;
671:   PetscScalar    *a;
672:   PetscReal      xir[2];
673:   PetscInt       Nf, p;
674:   const PetscInt dof = ctx->dof;

678:   DMGetNumFields(dm, &Nf);
679:   if (Nf) {DMGetField(dm, 0, NULL, (PetscObject *) &fem);}
680:   DMGetCoordinatesLocal(dm, &coordsLocal);
681:   DMGetCoordinateDM(dm, &dmCoord);
682:   SNESCreate(PETSC_COMM_SELF, &snes);
683:   SNESSetOptionsPrefix(snes, "quad_interp_");
684:   VecCreate(PETSC_COMM_SELF, &r);
685:   VecSetSizes(r, 2, 2);
686:   VecSetType(r,dm->vectype);
687:   VecDuplicate(r, &ref);
688:   VecDuplicate(r, &real);
689:   MatCreate(PETSC_COMM_SELF, &J);
690:   MatSetSizes(J, 2, 2, 2, 2);
691:   MatSetType(J, MATSEQDENSE);
692:   MatSetUp(J);
693:   SNESSetFunction(snes, r, QuadMap_Private, NULL);
694:   SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
695:   SNESGetKSP(snes, &ksp);
696:   KSPGetPC(ksp, &pc);
697:   PCSetType(pc, PCLU);
698:   SNESSetFromOptions(snes);

700:   VecGetArrayRead(ctx->coords, &coords);
701:   VecGetArray(v, &a);
702:   PetscFECreateTabulation(fem, 1, 1, xir, 0, &T);
703:   for (p = 0; p < ctx->n; ++p) {
704:     PetscScalar *x = NULL, *vertices = NULL;
705:     PetscScalar *xi;
706:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

708:     /* Can make this do all points at once */
709:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
710:     if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %d", coordSize, 4*2);
711:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
712:     SNESSetFunction(snes, NULL, NULL, (void*) vertices);
713:     SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
714:     VecGetArray(real, &xi);
715:     xi[0]  = coords[p*ctx->dim+0];
716:     xi[1]  = coords[p*ctx->dim+1];
717:     VecRestoreArray(real, &xi);
718:     SNESSolve(snes, real, ref);
719:     VecGetArray(ref, &xi);
720:     xir[0] = PetscRealPart(xi[0]);
721:     xir[1] = PetscRealPart(xi[1]);
722:     if (4*dof != xSize) {
723:       PetscInt d;

725:       xir[0] = 2.0*xir[0] - 1.0; xir[1] = 2.0*xir[1] - 1.0;
726:       PetscFEComputeTabulation(fem, 1, xir, 0, T);
727:       for (comp = 0; comp < dof; ++comp) {
728:         a[p*dof+comp] = 0.0;
729:         for (d = 0; d < xSize/dof; ++d) {
730:           a[p*dof+comp] += x[d*dof+comp]*T->T[0][d*dof+comp];
731:         }
732:       }
733:     } else {
734:       for (comp = 0; comp < dof; ++comp)
735:         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];
736:     }
737:     VecRestoreArray(ref, &xi);
738:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
739:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
740:   }
741:   PetscTabulationDestroy(&T);
742:   VecRestoreArray(v, &a);
743:   VecRestoreArrayRead(ctx->coords, &coords);

745:   SNESDestroy(&snes);
746:   VecDestroy(&r);
747:   VecDestroy(&ref);
748:   VecDestroy(&real);
749:   MatDestroy(&J);
750:   return(0);
751: }

753: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
754: {
755:   const PetscScalar *vertices = (const PetscScalar*) ctx;
756:   const PetscScalar x0        = vertices[0];
757:   const PetscScalar y0        = vertices[1];
758:   const PetscScalar z0        = vertices[2];
759:   const PetscScalar x1        = vertices[9];
760:   const PetscScalar y1        = vertices[10];
761:   const PetscScalar z1        = vertices[11];
762:   const PetscScalar x2        = vertices[6];
763:   const PetscScalar y2        = vertices[7];
764:   const PetscScalar z2        = vertices[8];
765:   const PetscScalar x3        = vertices[3];
766:   const PetscScalar y3        = vertices[4];
767:   const PetscScalar z3        = vertices[5];
768:   const PetscScalar x4        = vertices[12];
769:   const PetscScalar y4        = vertices[13];
770:   const PetscScalar z4        = vertices[14];
771:   const PetscScalar x5        = vertices[15];
772:   const PetscScalar y5        = vertices[16];
773:   const PetscScalar z5        = vertices[17];
774:   const PetscScalar x6        = vertices[18];
775:   const PetscScalar y6        = vertices[19];
776:   const PetscScalar z6        = vertices[20];
777:   const PetscScalar x7        = vertices[21];
778:   const PetscScalar y7        = vertices[22];
779:   const PetscScalar z7        = vertices[23];
780:   const PetscScalar f_1       = x1 - x0;
781:   const PetscScalar g_1       = y1 - y0;
782:   const PetscScalar h_1       = z1 - z0;
783:   const PetscScalar f_3       = x3 - x0;
784:   const PetscScalar g_3       = y3 - y0;
785:   const PetscScalar h_3       = z3 - z0;
786:   const PetscScalar f_4       = x4 - x0;
787:   const PetscScalar g_4       = y4 - y0;
788:   const PetscScalar h_4       = z4 - z0;
789:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
790:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
791:   const PetscScalar h_01      = z2 - z1 - z3 + z0;
792:   const PetscScalar f_12      = x7 - x3 - x4 + x0;
793:   const PetscScalar g_12      = y7 - y3 - y4 + y0;
794:   const PetscScalar h_12      = z7 - z3 - z4 + z0;
795:   const PetscScalar f_02      = x5 - x1 - x4 + x0;
796:   const PetscScalar g_02      = y5 - y1 - y4 + y0;
797:   const PetscScalar h_02      = z5 - z1 - z4 + z0;
798:   const PetscScalar f_012     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
799:   const PetscScalar g_012     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
800:   const PetscScalar h_012     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
801:   const PetscScalar *ref;
802:   PetscScalar       *real;
803:   PetscErrorCode    ierr;

806:   VecGetArrayRead(Xref,  &ref);
807:   VecGetArray(Xreal, &real);
808:   {
809:     const PetscScalar p0 = ref[0];
810:     const PetscScalar p1 = ref[1];
811:     const PetscScalar p2 = ref[2];

813:     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;
814:     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;
815:     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;
816:   }
817:   PetscLogFlops(114);
818:   VecRestoreArrayRead(Xref,  &ref);
819:   VecRestoreArray(Xreal, &real);
820:   return(0);
821: }

823: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
824: {
825:   const PetscScalar *vertices = (const PetscScalar*) ctx;
826:   const PetscScalar x0        = vertices[0];
827:   const PetscScalar y0        = vertices[1];
828:   const PetscScalar z0        = vertices[2];
829:   const PetscScalar x1        = vertices[9];
830:   const PetscScalar y1        = vertices[10];
831:   const PetscScalar z1        = vertices[11];
832:   const PetscScalar x2        = vertices[6];
833:   const PetscScalar y2        = vertices[7];
834:   const PetscScalar z2        = vertices[8];
835:   const PetscScalar x3        = vertices[3];
836:   const PetscScalar y3        = vertices[4];
837:   const PetscScalar z3        = vertices[5];
838:   const PetscScalar x4        = vertices[12];
839:   const PetscScalar y4        = vertices[13];
840:   const PetscScalar z4        = vertices[14];
841:   const PetscScalar x5        = vertices[15];
842:   const PetscScalar y5        = vertices[16];
843:   const PetscScalar z5        = vertices[17];
844:   const PetscScalar x6        = vertices[18];
845:   const PetscScalar y6        = vertices[19];
846:   const PetscScalar z6        = vertices[20];
847:   const PetscScalar x7        = vertices[21];
848:   const PetscScalar y7        = vertices[22];
849:   const PetscScalar z7        = vertices[23];
850:   const PetscScalar f_xy      = x2 - x1 - x3 + x0;
851:   const PetscScalar g_xy      = y2 - y1 - y3 + y0;
852:   const PetscScalar h_xy      = z2 - z1 - z3 + z0;
853:   const PetscScalar f_yz      = x7 - x3 - x4 + x0;
854:   const PetscScalar g_yz      = y7 - y3 - y4 + y0;
855:   const PetscScalar h_yz      = z7 - z3 - z4 + z0;
856:   const PetscScalar f_xz      = x5 - x1 - x4 + x0;
857:   const PetscScalar g_xz      = y5 - y1 - y4 + y0;
858:   const PetscScalar h_xz      = z5 - z1 - z4 + z0;
859:   const PetscScalar f_xyz     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
860:   const PetscScalar g_xyz     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
861:   const PetscScalar h_xyz     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
862:   const PetscScalar *ref;
863:   PetscErrorCode    ierr;

866:   VecGetArrayRead(Xref,  &ref);
867:   {
868:     const PetscScalar x       = ref[0];
869:     const PetscScalar y       = ref[1];
870:     const PetscScalar z       = ref[2];
871:     const PetscInt    rows[3] = {0, 1, 2};
872:     PetscScalar       values[9];

874:     values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
875:     values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
876:     values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
877:     values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
878:     values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
879:     values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
880:     values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
881:     values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
882:     values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;

884:     MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
885:   }
886:   PetscLogFlops(152);
887:   VecRestoreArrayRead(Xref,  &ref);
888:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
889:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
890:   return(0);
891: }

893: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
894: {
895:   DM             dmCoord;
896:   SNES           snes;
897:   KSP            ksp;
898:   PC             pc;
899:   Vec            coordsLocal, r, ref, real;
900:   Mat            J;
901:   const PetscScalar *coords;
902:   PetscScalar    *a;
903:   PetscInt       p;

907:   DMGetCoordinatesLocal(dm, &coordsLocal);
908:   DMGetCoordinateDM(dm, &dmCoord);
909:   SNESCreate(PETSC_COMM_SELF, &snes);
910:   SNESSetOptionsPrefix(snes, "hex_interp_");
911:   VecCreate(PETSC_COMM_SELF, &r);
912:   VecSetSizes(r, 3, 3);
913:   VecSetType(r,dm->vectype);
914:   VecDuplicate(r, &ref);
915:   VecDuplicate(r, &real);
916:   MatCreate(PETSC_COMM_SELF, &J);
917:   MatSetSizes(J, 3, 3, 3, 3);
918:   MatSetType(J, MATSEQDENSE);
919:   MatSetUp(J);
920:   SNESSetFunction(snes, r, HexMap_Private, NULL);
921:   SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
922:   SNESGetKSP(snes, &ksp);
923:   KSPGetPC(ksp, &pc);
924:   PCSetType(pc, PCLU);
925:   SNESSetFromOptions(snes);

927:   VecGetArrayRead(ctx->coords, &coords);
928:   VecGetArray(v, &a);
929:   for (p = 0; p < ctx->n; ++p) {
930:     PetscScalar *x = NULL, *vertices = NULL;
931:     PetscScalar *xi;
932:     PetscReal    xir[3];
933:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

935:     /* Can make this do all points at once */
936:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
937:     if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %d", coordSize, 8*3);
938:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
939:     if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %D", xSize, 8*ctx->dof);
940:     SNESSetFunction(snes, NULL, NULL, (void*) vertices);
941:     SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
942:     VecGetArray(real, &xi);
943:     xi[0]  = coords[p*ctx->dim+0];
944:     xi[1]  = coords[p*ctx->dim+1];
945:     xi[2]  = coords[p*ctx->dim+2];
946:     VecRestoreArray(real, &xi);
947:     SNESSolve(snes, real, ref);
948:     VecGetArray(ref, &xi);
949:     xir[0] = PetscRealPart(xi[0]);
950:     xir[1] = PetscRealPart(xi[1]);
951:     xir[2] = PetscRealPart(xi[2]);
952:     for (comp = 0; comp < ctx->dof; ++comp) {
953:       a[p*ctx->dof+comp] =
954:         x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
955:         x[3*ctx->dof+comp]*    xir[0]*(1-xir[1])*(1-xir[2]) +
956:         x[2*ctx->dof+comp]*    xir[0]*    xir[1]*(1-xir[2]) +
957:         x[1*ctx->dof+comp]*(1-xir[0])*    xir[1]*(1-xir[2]) +
958:         x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*   xir[2] +
959:         x[5*ctx->dof+comp]*    xir[0]*(1-xir[1])*   xir[2] +
960:         x[6*ctx->dof+comp]*    xir[0]*    xir[1]*   xir[2] +
961:         x[7*ctx->dof+comp]*(1-xir[0])*    xir[1]*   xir[2];
962:     }
963:     VecRestoreArray(ref, &xi);
964:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
965:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
966:   }
967:   VecRestoreArray(v, &a);
968:   VecRestoreArrayRead(ctx->coords, &coords);

970:   SNESDestroy(&snes);
971:   VecDestroy(&r);
972:   VecDestroy(&ref);
973:   VecDestroy(&real);
974:   MatDestroy(&J);
975:   return(0);
976: }

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

981:   Input Parameters:
982: + ctx - The DMInterpolationInfo context
983: . dm  - The DM
984: - x   - The local vector containing the field to be interpolated

986:   Output Parameters:
987: . v   - The vector containing the interpolated values

989:   Note: A suitable v can be obtained using DMInterpolationGetVector().

991:   Level: beginner

993: .seealso: DMInterpolationGetVector(), DMInterpolationAddPoints(), DMInterpolationCreate()
994: @*/
995: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
996: {
997:   PetscInt       dim, coneSize, n;

1004:   VecGetLocalSize(v, &n);
1005:   if (n != ctx->n*ctx->dof) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %D should be %D", n, ctx->n*ctx->dof);
1006:   if (n) {
1007:     DMGetDimension(dm, &dim);
1008:     DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);
1009:     if (dim == 2) {
1010:       if (coneSize == 3) {
1011:         DMInterpolate_Triangle_Private(ctx, dm, x, v);
1012:       } else if (coneSize == 4) {
1013:         DMInterpolate_Quad_Private(ctx, dm, x, v);
1014:       } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %D for point interpolation", dim);
1015:     } else if (dim == 3) {
1016:       if (coneSize == 4) {
1017:         DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
1018:       } else {
1019:         DMInterpolate_Hex_Private(ctx, dm, x, v);
1020:       }
1021:     } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %D for point interpolation", dim);
1022:   }
1023:   return(0);
1024: }

1026: /*@C
1027:   DMInterpolationDestroy - Destroys a DMInterpolationInfo context

1029:   Collective on ctx

1031:   Input Parameter:
1032: . ctx - the context

1034:   Level: beginner

1036: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
1037: @*/
1038: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1039: {

1044:   VecDestroy(&(*ctx)->coords);
1045:   PetscFree((*ctx)->points);
1046:   PetscFree((*ctx)->cells);
1047:   PetscFree(*ctx);
1048:   *ctx = NULL;
1049:   return(0);
1050: }

1052: /*@C
1053:   SNESMonitorFields - Monitors the residual for each field separately

1055:   Collective on SNES

1057:   Input Parameters:
1058: + snes   - the SNES context
1059: . its    - iteration number
1060: . fgnorm - 2-norm of residual
1061: - vf  - PetscViewerAndFormat of type ASCII

1063:   Notes:
1064:   This routine prints the residual norm at each iteration.

1066:   Level: intermediate

1068: .seealso: SNESMonitorSet(), SNESMonitorDefault()
1069: @*/
1070: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1071: {
1072:   PetscViewer        viewer = vf->viewer;
1073:   Vec                res;
1074:   DM                 dm;
1075:   PetscSection       s;
1076:   const PetscScalar *r;
1077:   PetscReal         *lnorms, *norms;
1078:   PetscInt           numFields, f, pStart, pEnd, p;
1079:   PetscErrorCode     ierr;

1083:   SNESGetFunction(snes, &res, NULL, NULL);
1084:   SNESGetDM(snes, &dm);
1085:   DMGetLocalSection(dm, &s);
1086:   PetscSectionGetNumFields(s, &numFields);
1087:   PetscSectionGetChart(s, &pStart, &pEnd);
1088:   PetscCalloc2(numFields, &lnorms, numFields, &norms);
1089:   VecGetArrayRead(res, &r);
1090:   for (p = pStart; p < pEnd; ++p) {
1091:     for (f = 0; f < numFields; ++f) {
1092:       PetscInt fdof, foff, d;

1094:       PetscSectionGetFieldDof(s, p, f, &fdof);
1095:       PetscSectionGetFieldOffset(s, p, f, &foff);
1096:       for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
1097:     }
1098:   }
1099:   VecRestoreArrayRead(res, &r);
1100:   MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1101:   PetscViewerPushFormat(viewer,vf->format);
1102:   PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
1103:   PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
1104:   for (f = 0; f < numFields; ++f) {
1105:     if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
1106:     PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
1107:   }
1108:   PetscViewerASCIIPrintf(viewer, "]\n");
1109:   PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
1110:   PetscViewerPopFormat(viewer);
1111:   PetscFree2(lnorms, norms);
1112:   return(0);
1113: }

1115: /********************* Residual Computation **************************/

1117: /*@
1118:   DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user

1120:   Input Parameters:
1121: + dm - The mesh
1122: . X  - Local solution
1123: - user - The user context

1125:   Output Parameter:
1126: . F  - Local output vector

1128:   Level: developer

1130: .seealso: DMPlexComputeJacobianAction()
1131: @*/
1132: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1133: {
1134:   DM             plex;
1135:   IS             allcellIS;
1136:   PetscInt       Nds, s, depth;

1140:   DMGetNumDS(dm, &Nds);
1141:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1142:   DMPlexGetDepth(plex, &depth);
1143:   DMGetStratumIS(plex, "dim", depth, &allcellIS);
1144:   if (!allcellIS) {DMGetStratumIS(plex, "depth", depth, &allcellIS);}
1145:   for (s = 0; s < Nds; ++s) {
1146:     PetscDS ds;
1147:     DMLabel label;
1148:     IS      cellIS;

1150:     DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1151:     if (!label) {
1152:       PetscObjectReference((PetscObject) allcellIS);
1153:       cellIS = allcellIS;
1154:     } else {
1155:       IS pointIS;

1157:       DMLabelGetStratumIS(label, 1, &pointIS);
1158:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1159:       ISDestroy(&pointIS);
1160:     }
1161:     DMPlexComputeResidual_Internal(plex, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1162:     ISDestroy(&cellIS);
1163:   }
1164:   ISDestroy(&allcellIS);
1165:   DMDestroy(&plex);
1166:   return(0);
1167: }

1169: /*@
1170:   DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X

1172:   Input Parameters:
1173: + dm - The mesh
1174: - user - The user context

1176:   Output Parameter:
1177: . X  - Local solution

1179:   Level: developer

1181: .seealso: DMPlexComputeJacobianAction()
1182: @*/
1183: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1184: {
1185:   DM             plex;

1189:   DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
1190:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
1191:   DMDestroy(&plex);
1192:   return(0);
1193: }

1195: /*@
1196:   DMPlexComputeJacobianAction - Form the local portion of the Jacobian action Z = J(X) Y at the local solution X using pointwise functions specified by the user.

1198:   Input Parameters:
1199: + dm - The mesh
1200: . cellIS -
1201: . t  - The time
1202: . X_tShift - The multiplier for the Jacobian with repsect to X_t
1203: . X  - Local solution vector
1204: . X_t  - Time-derivative of the local solution vector
1205: . Y  - Local input vector
1206: - user - The user context

1208:   Output Parameter:
1209: . Z - Local output vector

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

1215:   Level: developer

1217: .seealso: FormFunctionLocal()
1218: @*/
1219: PetscErrorCode DMPlexComputeJacobianAction(DM dm, IS cellIS, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Vec Y, Vec Z, void *user)
1220: {
1221:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1222:   const char       *name  = "Jacobian";
1223:   DM                dmAux, plex, plexAux = NULL;
1224:   DMEnclosureType   encAux;
1225:   Vec               A;
1226:   PetscDS           prob, probAux = NULL;
1227:   PetscQuadrature   quad;
1228:   PetscSection      section, globalSection, sectionAux;
1229:   PetscScalar      *elemMat, *elemMatD, *u, *u_t, *a = NULL, *y, *z;
1230:   PetscInt          Nf, fieldI, fieldJ;
1231:   PetscInt          totDim, totDimAux = 0;
1232:   const PetscInt   *cells;
1233:   PetscInt          cStart, cEnd, numCells, c;
1234:   PetscBool         hasDyn;
1235:   DMField           coordField;
1236:   PetscErrorCode    ierr;

1239:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
1240:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1241:   if (!cellIS) {
1242:     PetscInt depth;

1244:     DMPlexGetDepth(plex, &depth);
1245:     DMGetStratumIS(plex, "dim", depth, &cellIS);
1246:     if (!cellIS) {DMGetStratumIS(plex, "depth", depth, &cellIS);}
1247:   } else {
1248:     PetscObjectReference((PetscObject) cellIS);
1249:   }
1250:   DMGetLocalSection(dm, &section);
1251:   DMGetGlobalSection(dm, &globalSection);
1252:   ISGetLocalSize(cellIS, &numCells);
1253:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
1254:   DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);
1255:   PetscDSGetNumFields(prob, &Nf);
1256:   PetscDSGetTotalDimension(prob, &totDim);
1257:   PetscDSHasDynamicJacobian(prob, &hasDyn);
1258:   hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
1259:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1260:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1261:   if (dmAux) {
1262:     DMGetEnclosureRelation(dmAux, dm, &encAux);
1263:     DMConvert(dmAux, DMPLEX, &plexAux);
1264:     DMGetLocalSection(plexAux, &sectionAux);
1265:     DMGetDS(dmAux, &probAux);
1266:     PetscDSGetTotalDimension(probAux, &totDimAux);
1267:   }
1268:   VecSet(Z, 0.0);
1269:   PetscMalloc6(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat,hasDyn ? numCells*totDim*totDim : 0, &elemMatD,numCells*totDim,&y,totDim,&z);
1270:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1271:   DMGetCoordinateField(dm, &coordField);
1272:   for (c = cStart; c < cEnd; ++c) {
1273:     const PetscInt cell = cells ? cells[c] : c;
1274:     const PetscInt cind = c - cStart;
1275:     PetscScalar   *x = NULL,  *x_t = NULL;
1276:     PetscInt       i;

1278:     DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
1279:     for (i = 0; i < totDim; ++i) u[cind*totDim+i] = x[i];
1280:     DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
1281:     if (X_t) {
1282:       DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
1283:       for (i = 0; i < totDim; ++i) u_t[cind*totDim+i] = x_t[i];
1284:       DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
1285:     }
1286:     if (dmAux) {
1287:       PetscInt subcell;
1288:       DMGetEnclosurePoint(dmAux, dm, encAux, cell, &subcell);
1289:       DMPlexVecGetClosure(plexAux, sectionAux, A, subcell, NULL, &x);
1290:       for (i = 0; i < totDimAux; ++i) a[cind*totDimAux+i] = x[i];
1291:       DMPlexVecRestoreClosure(plexAux, sectionAux, A, subcell, NULL, &x);
1292:     }
1293:     DMPlexVecGetClosure(dm, section, Y, cell, NULL, &x);
1294:     for (i = 0; i < totDim; ++i) y[cind*totDim+i] = x[i];
1295:     DMPlexVecRestoreClosure(dm, section, Y, cell, NULL, &x);
1296:   }
1297:   PetscArrayzero(elemMat, numCells*totDim*totDim);
1298:   if (hasDyn)  {PetscArrayzero(elemMatD, numCells*totDim*totDim);}
1299:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1300:     PetscFE  fe;
1301:     PetscInt Nb;
1302:     /* Conforming batches */
1303:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1304:     /* Remainder */
1305:     PetscInt Nr, offset, Nq;
1306:     PetscQuadrature qGeom = NULL;
1307:     PetscInt    maxDegree;
1308:     PetscFEGeom *cgeomFEM, *chunkGeom = NULL, *remGeom = NULL;

1310:     PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1311:     PetscFEGetQuadrature(fe, &quad);
1312:     PetscFEGetDimension(fe, &Nb);
1313:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1314:     DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1315:     if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField,cellIS,&qGeom);}
1316:     if (!qGeom) {
1317:       PetscFEGetQuadrature(fe,&qGeom);
1318:       PetscObjectReference((PetscObject)qGeom);
1319:     }
1320:     PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
1321:     DMSNESGetFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);
1322:     blockSize = Nb;
1323:     batchSize = numBlocks * blockSize;
1324:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1325:     numChunks = numCells / (numBatches*batchSize);
1326:     Ne        = numChunks*numBatches*batchSize;
1327:     Nr        = numCells % (numBatches*batchSize);
1328:     offset    = numCells - Nr;
1329:     PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1330:     PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&remGeom);
1331:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1332:       PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMat);
1333:       PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMat[offset*totDim*totDim]);
1334:       if (hasDyn) {
1335:         PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMatD);
1336:         PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMatD[offset*totDim*totDim]);
1337:       }
1338:     }
1339:     PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&remGeom);
1340:     PetscFEGeomRestoreChunk(cgeomFEM,0,offset,&chunkGeom);
1341:     DMSNESRestoreFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);
1342:     PetscQuadratureDestroy(&qGeom);
1343:   }
1344:   if (hasDyn) {
1345:     for (c = 0; c < numCells*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
1346:   }
1347:   for (c = cStart; c < cEnd; ++c) {
1348:     const PetscInt     cell = cells ? cells[c] : c;
1349:     const PetscInt     cind = c - cStart;
1350:     const PetscBLASInt M = totDim, one = 1;
1351:     const PetscScalar  a = 1.0, b = 0.0;

1353:     PetscStackCallBLAS("BLASgemv", BLASgemv_("N", &M, &M, &a, &elemMat[cind*totDim*totDim], &M, &y[cind*totDim], &one, &b, z, &one));
1354:     if (mesh->printFEM > 1) {
1355:       DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[cind*totDim*totDim]);
1356:       DMPrintCellVector(c, "Y",  totDim, &y[cind*totDim]);
1357:       DMPrintCellVector(c, "Z",  totDim, z);
1358:     }
1359:     DMPlexVecSetClosure(dm, section, Z, cell, z, ADD_VALUES);
1360:   }
1361:   PetscFree6(u,u_t,elemMat,elemMatD,y,z);
1362:   if (mesh->printFEM) {
1363:     PetscPrintf(PetscObjectComm((PetscObject)Z), "Z:\n");
1364:     VecView(Z, NULL);
1365:   }
1366:   PetscFree(a);
1367:   ISDestroy(&cellIS);
1368:   DMDestroy(&plexAux);
1369:   DMDestroy(&plex);
1370:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
1371:   return(0);
1372: }

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

1377:   Input Parameters:
1378: + dm - The mesh
1379: . X  - Local input vector
1380: - user - The user context

1382:   Output Parameter:
1383: . Jac  - Jacobian matrix

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

1389:   Level: developer

1391: .seealso: FormFunctionLocal()
1392: @*/
1393: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1394: {
1395:   DM             plex;
1396:   IS             allcellIS;
1397:   PetscBool      hasJac, hasPrec;
1398:   PetscInt       Nds, s, depth;

1402:   DMGetNumDS(dm, &Nds);
1403:   DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1404:   DMPlexGetDepth(plex, &depth);
1405:   DMGetStratumIS(plex, "dim", depth, &allcellIS);
1406:   if (!allcellIS) {DMGetStratumIS(plex, "depth", depth, &allcellIS);}
1407:   for (s = 0; s < Nds; ++s) {
1408:     PetscDS ds;
1409:     DMLabel label;
1410:     IS      cellIS;

1412:     DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1413:     if (!label) {
1414:       PetscObjectReference((PetscObject) allcellIS);
1415:       cellIS = allcellIS;
1416:     } else {
1417:       IS pointIS;

1419:       DMLabelGetStratumIS(label, 1, &pointIS);
1420:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1421:       ISDestroy(&pointIS);
1422:     }
1423:     if (!s) {
1424:       PetscDSHasJacobian(ds, &hasJac);
1425:       PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1426:       if (hasJac && hasPrec) {MatZeroEntries(Jac);}
1427:       MatZeroEntries(JacP);
1428:     }
1429:     DMPlexComputeJacobian_Internal(plex, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user);
1430:     ISDestroy(&cellIS);
1431:   }
1432:   ISDestroy(&allcellIS);
1433:   DMDestroy(&plex);
1434:   return(0);
1435: }

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

1440:    Input Parameters:
1441: +     X - SNES linearization point
1442: .     ovl - index set of overlapping subdomains

1444:    Output Parameter:
1445: .     J - unassembled (Neumann) local matrix

1447:    Level: intermediate

1449: .seealso: DMCreateNeumannOverlap(), MATIS, PCHPDDMSetAuxiliaryMat()
1450: */
1451: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1452: {
1453:   SNES           snes;
1454:   Mat            pJ;
1455:   DM             ovldm,origdm;
1456:   DMSNES         sdm;
1457:   PetscErrorCode (*bfun)(DM,Vec,void*);
1458:   PetscErrorCode (*jfun)(DM,Vec,Mat,Mat,void*);
1459:   void           *bctx,*jctx;

1463:   PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_MATIS",(PetscObject*)&pJ);
1464:   if (!pJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing overlapping Mat");
1465:   PetscObjectQuery((PetscObject)ovl,"_DM_Original_HPDDM",(PetscObject*)&origdm);
1466:   if (!origdm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing original DM");
1467:   MatGetDM(pJ,&ovldm);
1468:   DMSNESGetBoundaryLocal(origdm,&bfun,&bctx);
1469:   DMSNESSetBoundaryLocal(ovldm,bfun,bctx);
1470:   DMSNESGetJacobianLocal(origdm,&jfun,&jctx);
1471:   DMSNESSetJacobianLocal(ovldm,jfun,jctx);
1472:   PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject*)&snes);
1473:   if (!snes) {
1474:     SNESCreate(PetscObjectComm((PetscObject)ovl),&snes);
1475:     SNESSetDM(snes,ovldm);
1476:     PetscObjectCompose((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject)snes);
1477:     PetscObjectDereference((PetscObject)snes);
1478:   }
1479:   DMGetDMSNES(ovldm,&sdm);
1480:   VecLockReadPush(X);
1481:   PetscStackPush("SNES user Jacobian function");
1482:   (*sdm->ops->computejacobian)(snes,X,pJ,pJ,sdm->jacobianctx);
1483:   PetscStackPop;
1484:   VecLockReadPop(X);
1485:   /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1486:   {
1487:     Mat locpJ;

1489:     MatISGetLocalMat(pJ,&locpJ);
1490:     MatCopy(locpJ,J,SAME_NONZERO_PATTERN);
1491:   }
1492:   return(0);
1493: }

1495: /*@
1496:   DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian.

1498:   Input Parameters:
1499: + dm - The DM object
1500: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see PetscDSAddBoundary())
1501: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
1502: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())

1504:   Level: developer
1505: @*/
1506: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1507: {

1511:   DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);
1512:   DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);
1513:   DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);
1514:   PetscObjectComposeFunction((PetscObject)dm,"MatComputeNeumannOverlap_C",MatComputeNeumannOverlap_Plex);
1515:   return(0);
1516: }

1518: /*@C
1519:   DMSNESCheckDiscretization - Check the discretization error of the exact solution

1521:   Input Parameters:
1522: + snes - the SNES object
1523: . dm   - the DM
1524: . t    - the time
1525: . u    - a DM vector
1526: - tol  - A tolerance for the check, or -1 to print the results instead

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

1531:   Note: The user must call PetscDSSetExactSolution() beforehand

1533:   Level: developer

1535: .seealso: DNSNESCheckFromOptions(), DMSNESCheckResidual(), DMSNESCheckJacobian(), PetscDSSetExactSolution()
1536: @*/
1537: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1538: {
1539:   PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1540:   void            **ectxs;
1541:   PetscReal        *err;
1542:   MPI_Comm          comm;
1543:   PetscInt          Nf, f;
1544:   PetscErrorCode    ierr;


1552:   DMComputeExactSolution(dm, t, u, NULL);
1553:   VecViewFromOptions(u, NULL, "-vec_view");

1555:   PetscObjectGetComm((PetscObject) snes, &comm);
1556:   DMGetNumFields(dm, &Nf);
1557:   PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);
1558:   {
1559:     PetscInt Nds, s;

1561:     DMGetNumDS(dm, &Nds);
1562:     for (s = 0; s < Nds; ++s) {
1563:       PetscDS         ds;
1564:       DMLabel         label;
1565:       IS              fieldIS;
1566:       const PetscInt *fields;
1567:       PetscInt        dsNf, f;

1569:       DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);
1570:       PetscDSGetNumFields(ds, &dsNf);
1571:       ISGetIndices(fieldIS, &fields);
1572:       for (f = 0; f < dsNf; ++f) {
1573:         const PetscInt field = fields[f];
1574:         PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);
1575:       }
1576:       ISRestoreIndices(fieldIS, &fields);
1577:     }
1578:   }
1579:   if (Nf > 1) {
1580:     DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err);
1581:     if (tol >= 0.0) {
1582:       for (f = 0; f < Nf; ++f) {
1583:         if (err[f] > tol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g for field %D exceeds tolerance %g", (double) err[f], f, (double) tol);
1584:       }
1585:     } else if (error) {
1586:       for (f = 0; f < Nf; ++f) error[f] = err[f];
1587:     } else {
1588:       PetscPrintf(comm, "L_2 Error: [");
1589:       for (f = 0; f < Nf; ++f) {
1590:         if (f) {PetscPrintf(comm, ", ");}
1591:         PetscPrintf(comm, "%g", (double)err[f]);
1592:       }
1593:       PetscPrintf(comm, "]\n");
1594:     }
1595:   } else {
1596:     DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]);
1597:     if (tol >= 0.0) {
1598:       if (err[0] > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double) err[0], (double) tol);
1599:     } else if (error) {
1600:       error[0] = err[0];
1601:     } else {
1602:       PetscPrintf(comm, "L_2 Error: %g\n", (double) err[0]);
1603:     }
1604:   }
1605:   PetscFree3(exacts, ectxs, err);
1606:   return(0);
1607: }

1609: /*@C
1610:   DMSNESCheckResidual - Check the residual of the exact solution

1612:   Input Parameters:
1613: + snes - the SNES object
1614: . dm   - the DM
1615: . u    - a DM vector
1616: - tol  - A tolerance for the check, or -1 to print the results instead

1618:   Output Parameters:
1619: . residual - The residual norm of the exact solution, or NULL

1621:   Level: developer

1623: .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
1624: @*/
1625: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1626: {
1627:   MPI_Comm       comm;
1628:   Vec            r;
1629:   PetscReal      res;

1637:   PetscObjectGetComm((PetscObject) snes, &comm);
1638:   DMComputeExactSolution(dm, 0.0, u, NULL);
1639:   VecDuplicate(u, &r);
1640:   SNESComputeFunction(snes, u, r);
1641:   VecNorm(r, NORM_2, &res);
1642:   if (tol >= 0.0) {
1643:     if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
1644:   } else if (residual) {
1645:     *residual = res;
1646:   } else {
1647:     PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
1648:     VecChop(r, 1.0e-10);
1649:     PetscObjectSetName((PetscObject) r, "Initial Residual");
1650:     PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
1651:     VecViewFromOptions(r, NULL, "-vec_view");
1652:   }
1653:   VecDestroy(&r);
1654:   return(0);
1655: }

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

1660:   Input Parameters:
1661: + snes - the SNES object
1662: . dm   - the DM
1663: . u    - a DM vector
1664: - tol  - A tolerance for the check, or -1 to print the results instead

1666:   Output Parameters:
1667: + isLinear - Flag indicaing that the function looks linear, or NULL
1668: - convRate - The rate of convergence of the linear model, or NULL

1670:   Level: developer

1672: .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
1673: @*/
1674: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1675: {
1676:   MPI_Comm       comm;
1677:   PetscDS        ds;
1678:   Mat            J, M;
1679:   MatNullSpace   nullspace;
1680:   PetscReal      slope, intercept;
1681:   PetscBool      hasJac, hasPrec, isLin = PETSC_FALSE;

1690:   PetscObjectGetComm((PetscObject) snes, &comm);
1691:   DMComputeExactSolution(dm, 0.0, u, NULL);
1692:   /* Create and view matrices */
1693:   DMCreateMatrix(dm, &J);
1694:   DMGetDS(dm, &ds);
1695:   PetscDSHasJacobian(ds, &hasJac);
1696:   PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1697:   if (hasJac && hasPrec) {
1698:     DMCreateMatrix(dm, &M);
1699:     SNESComputeJacobian(snes, u, J, M);
1700:     PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
1701:     PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
1702:     MatViewFromOptions(M, NULL, "-mat_view");
1703:     MatDestroy(&M);
1704:   } else {
1705:     SNESComputeJacobian(snes, u, J, J);
1706:   }
1707:   PetscObjectSetName((PetscObject) J, "Jacobian");
1708:   PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
1709:   MatViewFromOptions(J, NULL, "-mat_view");
1710:   /* Check nullspace */
1711:   MatGetNullSpace(J, &nullspace);
1712:   if (nullspace) {
1713:     PetscBool isNull;
1714:     MatNullSpaceTest(nullspace, J, &isNull);
1715:     if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1716:   }
1717:   /* Taylor test */
1718:   {
1719:     PetscRandom rand;
1720:     Vec         du, uhat, r, rhat, df;
1721:     PetscReal   h;
1722:     PetscReal  *es, *hs, *errors;
1723:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1724:     PetscInt    Nv, v;

1726:     /* Choose a perturbation direction */
1727:     PetscRandomCreate(comm, &rand);
1728:     VecDuplicate(u, &du);
1729:     VecSetRandom(du, rand);
1730:     PetscRandomDestroy(&rand);
1731:     VecDuplicate(u, &df);
1732:     MatMult(J, du, df);
1733:     /* Evaluate residual at u, F(u), save in vector r */
1734:     VecDuplicate(u, &r);
1735:     SNESComputeFunction(snes, u, r);
1736:     /* Look at the convergence of our Taylor approximation as we approach u */
1737:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1738:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
1739:     VecDuplicate(u, &uhat);
1740:     VecDuplicate(u, &rhat);
1741:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1742:       VecWAXPY(uhat, h, du, u);
1743:       /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1744:       SNESComputeFunction(snes, uhat, rhat);
1745:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
1746:       VecNorm(rhat, NORM_2, &errors[Nv]);

1748:       es[Nv] = PetscLog10Real(errors[Nv]);
1749:       hs[Nv] = PetscLog10Real(h);
1750:     }
1751:     VecDestroy(&uhat);
1752:     VecDestroy(&rhat);
1753:     VecDestroy(&df);
1754:     VecDestroy(&r);
1755:     VecDestroy(&du);
1756:     for (v = 0; v < Nv; ++v) {
1757:       if ((tol >= 0) && (errors[v] > tol)) break;
1758:       else if (errors[v] > PETSC_SMALL)    break;
1759:     }
1760:     if (v == Nv) isLin = PETSC_TRUE;
1761:     PetscLinearRegression(Nv, hs, es, &slope, &intercept);
1762:     PetscFree3(es, hs, errors);
1763:     /* Slope should be about 2 */
1764:     if (tol >= 0) {
1765:       if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
1766:     } else if (isLinear || convRate) {
1767:       if (isLinear) *isLinear = isLin;
1768:       if (convRate) *convRate = slope;
1769:     } else {
1770:       if (!isLin) {PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);}
1771:       else        {PetscPrintf(comm, "Function appears to be linear\n");}
1772:     }
1773:   }
1774:   MatDestroy(&J);
1775:   return(0);
1776: }

1778: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1779: {

1783:   DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL);
1784:   DMSNESCheckResidual(snes, dm, u, -1.0, NULL);
1785:   DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL);
1786:   return(0);
1787: }

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

1792:   Input Parameters:
1793: + snes - the SNES object
1794: - u    - representative SNES vector

1796:   Note: The user must call PetscDSSetExactSolution() beforehand

1798:   Level: developer
1799: @*/
1800: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1801: {
1802:   DM             dm;
1803:   Vec            sol;
1804:   PetscBool      check;

1808:   PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);
1809:   if (!check) return(0);
1810:   SNESGetDM(snes, &dm);
1811:   VecDuplicate(u, &sol);
1812:   SNESSetSolution(snes, sol);
1813:   DMSNESCheck_Internal(snes, dm, sol);
1814:   VecDestroy(&sol);
1815:   return(0);
1816: }