Actual source code: dmplexsnes.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I "petscdmplex.h" I*/
  2: #include <petsc/private/snesimpl.h>     /*I "petscsnes.h"   I*/
  3: #include <petscds.h>
  4: #include <petsc/private/petscimpl.h>
  5: #include <petsc/private/petscfeimpl.h>

  7: /************************** Interpolation *******************************/

 11: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
 12: {
 13:   PetscBool      isPlex;

 17:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
 18:   if (isPlex) {
 19:     *plex = dm;
 20:     PetscObjectReference((PetscObject) dm);
 21:   } else {
 22:     PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
 23:     if (!*plex) {
 24:       DMConvert(dm,DMPLEX,plex);
 25:       PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
 26:       if (copy) {
 27:         PetscInt    i;
 28:         PetscObject obj;
 29:         const char *comps[3] = {"A","dmAux","dmCh"};

 31:         DMCopyDMSNES(dm, *plex);
 32:         for (i = 0; i < 3; i++) {
 33:           PetscObjectQuery((PetscObject) dm, comps[i], &obj);
 34:           PetscObjectCompose((PetscObject) *plex, comps[i], obj);
 35:         }
 36:       }
 37:     } else {
 38:       PetscObjectReference((PetscObject) *plex);
 39:     }
 40:   }
 41:   return(0);
 42: }

 46: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
 47: {

 52:   PetscMalloc(sizeof(struct _DMInterpolationInfo), ctx);

 54:   (*ctx)->comm   = comm;
 55:   (*ctx)->dim    = -1;
 56:   (*ctx)->nInput = 0;
 57:   (*ctx)->points = NULL;
 58:   (*ctx)->cells  = NULL;
 59:   (*ctx)->n      = -1;
 60:   (*ctx)->coords = NULL;
 61:   return(0);
 62: }

 66: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
 67: {
 69:   if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);
 70:   ctx->dim = dim;
 71:   return(0);
 72: }

 76: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
 77: {
 80:   *dim = ctx->dim;
 81:   return(0);
 82: }

 86: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
 87: {
 89:   if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);
 90:   ctx->dof = dof;
 91:   return(0);
 92: }

 96: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
 97: {
100:   *dof = ctx->dof;
101:   return(0);
102: }

106: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
107: {

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

115:   PetscMalloc1(n*ctx->dim, &ctx->points);
116:   PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));
117:   return(0);
118: }

122: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
123: {
124:   MPI_Comm          comm = ctx->comm;
125:   PetscScalar       *a;
126:   PetscInt          p, q, i;
127:   PetscMPIInt       rank, size;
128:   PetscErrorCode    ierr;
129:   Vec               pointVec;
130:   PetscSF           cellSF;
131:   PetscLayout       layout;
132:   PetscReal         *globalPoints;
133:   PetscScalar       *globalPointsScalar;
134:   const PetscInt    *ranges;
135:   PetscMPIInt       *counts, *displs;
136:   const PetscSFNode *foundCells;
137:   const PetscInt    *foundPoints;
138:   PetscMPIInt       *foundProcs, *globalProcs;
139:   PetscInt          n, N, numFound;

143:   MPI_Comm_size(comm, &size);
144:   MPI_Comm_rank(comm, &rank);
145:   if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
146:   /* Locate points */
147:   n = ctx->nInput;
148:   if (!redundantPoints) {
149:     PetscLayoutCreate(comm, &layout);
150:     PetscLayoutSetBlockSize(layout, 1);
151:     PetscLayoutSetLocalSize(layout, n);
152:     PetscLayoutSetUp(layout);
153:     PetscLayoutGetSize(layout, &N);
154:     /* Communicate all points to all processes */
155:     PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
156:     PetscLayoutGetRanges(layout, &ranges);
157:     for (p = 0; p < size; ++p) {
158:       counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
159:       displs[p] = ranges[p]*ctx->dim;
160:     }
161:     MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
162:   } else {
163:     N = n;
164:     globalPoints = ctx->points;
165:     counts = displs = NULL;
166:     layout = NULL;
167:   }
168: #if 0
169:   PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
170:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
171: #else
172: #if defined(PETSC_USE_COMPLEX)
173:   PetscMalloc1(N,&globalPointsScalar);
174:   for (i=0; i<N; i++) globalPointsScalar[i] = globalPoints[i];
175: #else
176:   globalPointsScalar = globalPoints;
177: #endif
178:   VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
179:   PetscMalloc2(N,&foundProcs,N,&globalProcs);
180:   cellSF = NULL;
181:   DMLocatePoints(dm, pointVec, &cellSF);
182:   PetscSFGetGraph(cellSF,NULL,&numFound,&foundPoints,&foundCells);
183: #endif
184:   for (p = 0; p < numFound; ++p) {
185:     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
186:     else foundProcs[foundPoints ? foundPoints[p] : p] = size;
187:   }
188:   /* Let the lowest rank process own each point */
189:   MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
190:   ctx->n = 0;
191:   for (p = 0; p < N; ++p) {
192:     if (globalProcs[p] == size) SETERRQ4(comm, PETSC_ERR_PLIB, "Point %d: %g %g %g not located in mesh", p, globalPoints[p*ctx->dim+0], ctx->dim > 1 ? globalPoints[p*ctx->dim+1] : 0.0, ctx->dim > 2 ? globalPoints[p*ctx->dim+2] : 0.0);
193:     else if (globalProcs[p] == rank) ctx->n++;
194:   }
195:   /* Create coordinates vector and array of owned cells */
196:   PetscMalloc1(ctx->n, &ctx->cells);
197:   VecCreate(comm, &ctx->coords);
198:   VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
199:   VecSetBlockSize(ctx->coords, ctx->dim);
200:   VecSetType(ctx->coords,VECSTANDARD);
201:   VecGetArray(ctx->coords, &a);
202:   for (p = 0, q = 0, i = 0; p < N; ++p) {
203:     if (globalProcs[p] == rank) {
204:       PetscInt d;

206:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
207:       ctx->cells[q++] = foundCells[p].index;
208:     }
209:   }
210:   VecRestoreArray(ctx->coords, &a);
211: #if 0
212:   PetscFree3(foundCells,foundProcs,globalProcs);
213: #else
214:   PetscFree2(foundProcs,globalProcs);
215:   PetscSFDestroy(&cellSF);
216:   VecDestroy(&pointVec);
217: #endif
218:   if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
219:   if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
220:   PetscLayoutDestroy(&layout);
221:   return(0);
222: }

226: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
227: {
230:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
231:   *coordinates = ctx->coords;
232:   return(0);
233: }

237: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
238: {

243:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
244:   VecCreate(ctx->comm, v);
245:   VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
246:   VecSetBlockSize(*v, ctx->dof);
247:   VecSetType(*v,VECSTANDARD);
248:   return(0);
249: }

253: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
254: {

259:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
260:   VecDestroy(v);
261:   return(0);
262: }

266: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
267: {
268:   PetscReal      *v0, *J, *invJ, detJ;
269:   const PetscScalar *coords;
270:   PetscScalar    *a;
271:   PetscInt       p;

275:   PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
276:   VecGetArrayRead(ctx->coords, &coords);
277:   VecGetArray(v, &a);
278:   for (p = 0; p < ctx->n; ++p) {
279:     PetscInt     c = ctx->cells[p];
280:     PetscScalar *x = NULL;
281:     PetscReal    xi[4];
282:     PetscInt     d, f, comp;

284:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
285:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
286:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
287:     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];

289:     for (d = 0; d < ctx->dim; ++d) {
290:       xi[d] = 0.0;
291:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
292:       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];
293:     }
294:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
295:   }
296:   VecRestoreArray(v, &a);
297:   VecRestoreArrayRead(ctx->coords, &coords);
298:   PetscFree3(v0, J, invJ);
299:   return(0);
300: }

304: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
305: {
306:   PetscReal      *v0, *J, *invJ, detJ;
307:   const PetscScalar *coords;
308:   PetscScalar    *a;
309:   PetscInt       p;

313:   PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
314:   VecGetArrayRead(ctx->coords, &coords);
315:   VecGetArray(v, &a);
316:   for (p = 0; p < ctx->n; ++p) {
317:     PetscInt       c = ctx->cells[p];
318:     const PetscInt order[3] = {2, 1, 3};
319:     PetscScalar   *x = NULL;
320:     PetscReal      xi[4];
321:     PetscInt       d, f, comp;

323:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
324:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
325:     DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
326:     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];

328:     for (d = 0; d < ctx->dim; ++d) {
329:       xi[d] = 0.0;
330:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
331:       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];
332:     }
333:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
334:   }
335:   VecRestoreArray(v, &a);
336:   VecRestoreArrayRead(ctx->coords, &coords);
337:   PetscFree3(v0, J, invJ);
338:   return(0);
339: }

343: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
344: {
345:   const PetscScalar *vertices = (const PetscScalar*) ctx;
346:   const PetscScalar x0        = vertices[0];
347:   const PetscScalar y0        = vertices[1];
348:   const PetscScalar x1        = vertices[2];
349:   const PetscScalar y1        = vertices[3];
350:   const PetscScalar x2        = vertices[4];
351:   const PetscScalar y2        = vertices[5];
352:   const PetscScalar x3        = vertices[6];
353:   const PetscScalar y3        = vertices[7];
354:   const PetscScalar f_1       = x1 - x0;
355:   const PetscScalar g_1       = y1 - y0;
356:   const PetscScalar f_3       = x3 - x0;
357:   const PetscScalar g_3       = y3 - y0;
358:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
359:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
360:   const PetscScalar *ref;
361:   PetscScalar       *real;
362:   PetscErrorCode    ierr;

365:   VecGetArrayRead(Xref,  &ref);
366:   VecGetArray(Xreal, &real);
367:   {
368:     const PetscScalar p0 = ref[0];
369:     const PetscScalar p1 = ref[1];

371:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
372:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
373:   }
374:   PetscLogFlops(28);
375:   VecRestoreArrayRead(Xref,  &ref);
376:   VecRestoreArray(Xreal, &real);
377:   return(0);
378: }

380: #include <petsc/private/dmimpl.h>
383: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
384: {
385:   const PetscScalar *vertices = (const PetscScalar*) ctx;
386:   const PetscScalar x0        = vertices[0];
387:   const PetscScalar y0        = vertices[1];
388:   const PetscScalar x1        = vertices[2];
389:   const PetscScalar y1        = vertices[3];
390:   const PetscScalar x2        = vertices[4];
391:   const PetscScalar y2        = vertices[5];
392:   const PetscScalar x3        = vertices[6];
393:   const PetscScalar y3        = vertices[7];
394:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
395:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
396:   const PetscScalar *ref;
397:   PetscErrorCode    ierr;

400:   VecGetArrayRead(Xref,  &ref);
401:   {
402:     const PetscScalar x       = ref[0];
403:     const PetscScalar y       = ref[1];
404:     const PetscInt    rows[2] = {0, 1};
405:     PetscScalar       values[4];

407:     values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
408:     values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
409:     MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
410:   }
411:   PetscLogFlops(30);
412:   VecRestoreArrayRead(Xref,  &ref);
413:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
414:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
415:   return(0);
416: }

420: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
421: {
422:   DM             dmCoord;
423:   SNES           snes;
424:   KSP            ksp;
425:   PC             pc;
426:   Vec            coordsLocal, r, ref, real;
427:   Mat            J;
428:   const PetscScalar *coords;
429:   PetscScalar    *a;
430:   PetscInt       p;

434:   DMGetCoordinatesLocal(dm, &coordsLocal);
435:   DMGetCoordinateDM(dm, &dmCoord);
436:   SNESCreate(PETSC_COMM_SELF, &snes);
437:   SNESSetOptionsPrefix(snes, "quad_interp_");
438:   VecCreate(PETSC_COMM_SELF, &r);
439:   VecSetSizes(r, 2, 2);
440:   VecSetType(r,dm->vectype);
441:   VecDuplicate(r, &ref);
442:   VecDuplicate(r, &real);
443:   MatCreate(PETSC_COMM_SELF, &J);
444:   MatSetSizes(J, 2, 2, 2, 2);
445:   MatSetType(J, MATSEQDENSE);
446:   MatSetUp(J);
447:   SNESSetFunction(snes, r, QuadMap_Private, NULL);
448:   SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
449:   SNESGetKSP(snes, &ksp);
450:   KSPGetPC(ksp, &pc);
451:   PCSetType(pc, PCLU);
452:   SNESSetFromOptions(snes);

454:   VecGetArrayRead(ctx->coords, &coords);
455:   VecGetArray(v, &a);
456:   for (p = 0; p < ctx->n; ++p) {
457:     PetscScalar *x = NULL, *vertices = NULL;
458:     PetscScalar *xi;
459:     PetscReal    xir[2];
460:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

462:     /* Can make this do all points at once */
463:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
464:     if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);
465:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
466:     if (4*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 4*ctx->dof);
467:     SNESSetFunction(snes, NULL, NULL, (void*) vertices);
468:     SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
469:     VecGetArray(real, &xi);
470:     xi[0]  = coords[p*ctx->dim+0];
471:     xi[1]  = coords[p*ctx->dim+1];
472:     VecRestoreArray(real, &xi);
473:     SNESSolve(snes, real, ref);
474:     VecGetArray(ref, &xi);
475:     xir[0] = PetscRealPart(xi[0]);
476:     xir[1] = PetscRealPart(xi[1]);
477:     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp]*(1 - xir[0])*(1 - xir[1]) + x[1*ctx->dof+comp]*xir[0]*(1 - xir[1]) + x[2*ctx->dof+comp]*xir[0]*xir[1] + x[3*ctx->dof+comp]*(1 - xir[0])*xir[1];

479:     VecRestoreArray(ref, &xi);
480:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
481:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
482:   }
483:   VecRestoreArray(v, &a);
484:   VecRestoreArrayRead(ctx->coords, &coords);

486:   SNESDestroy(&snes);
487:   VecDestroy(&r);
488:   VecDestroy(&ref);
489:   VecDestroy(&real);
490:   MatDestroy(&J);
491:   return(0);
492: }

496: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
497: {
498:   const PetscScalar *vertices = (const PetscScalar*) ctx;
499:   const PetscScalar x0        = vertices[0];
500:   const PetscScalar y0        = vertices[1];
501:   const PetscScalar z0        = vertices[2];
502:   const PetscScalar x1        = vertices[9];
503:   const PetscScalar y1        = vertices[10];
504:   const PetscScalar z1        = vertices[11];
505:   const PetscScalar x2        = vertices[6];
506:   const PetscScalar y2        = vertices[7];
507:   const PetscScalar z2        = vertices[8];
508:   const PetscScalar x3        = vertices[3];
509:   const PetscScalar y3        = vertices[4];
510:   const PetscScalar z3        = vertices[5];
511:   const PetscScalar x4        = vertices[12];
512:   const PetscScalar y4        = vertices[13];
513:   const PetscScalar z4        = vertices[14];
514:   const PetscScalar x5        = vertices[15];
515:   const PetscScalar y5        = vertices[16];
516:   const PetscScalar z5        = vertices[17];
517:   const PetscScalar x6        = vertices[18];
518:   const PetscScalar y6        = vertices[19];
519:   const PetscScalar z6        = vertices[20];
520:   const PetscScalar x7        = vertices[21];
521:   const PetscScalar y7        = vertices[22];
522:   const PetscScalar z7        = vertices[23];
523:   const PetscScalar f_1       = x1 - x0;
524:   const PetscScalar g_1       = y1 - y0;
525:   const PetscScalar h_1       = z1 - z0;
526:   const PetscScalar f_3       = x3 - x0;
527:   const PetscScalar g_3       = y3 - y0;
528:   const PetscScalar h_3       = z3 - z0;
529:   const PetscScalar f_4       = x4 - x0;
530:   const PetscScalar g_4       = y4 - y0;
531:   const PetscScalar h_4       = z4 - z0;
532:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
533:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
534:   const PetscScalar h_01      = z2 - z1 - z3 + z0;
535:   const PetscScalar f_12      = x7 - x3 - x4 + x0;
536:   const PetscScalar g_12      = y7 - y3 - y4 + y0;
537:   const PetscScalar h_12      = z7 - z3 - z4 + z0;
538:   const PetscScalar f_02      = x5 - x1 - x4 + x0;
539:   const PetscScalar g_02      = y5 - y1 - y4 + y0;
540:   const PetscScalar h_02      = z5 - z1 - z4 + z0;
541:   const PetscScalar f_012     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
542:   const PetscScalar g_012     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
543:   const PetscScalar h_012     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
544:   const PetscScalar *ref;
545:   PetscScalar       *real;
546:   PetscErrorCode    ierr;

549:   VecGetArrayRead(Xref,  &ref);
550:   VecGetArray(Xreal, &real);
551:   {
552:     const PetscScalar p0 = ref[0];
553:     const PetscScalar p1 = ref[1];
554:     const PetscScalar p2 = ref[2];

556:     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;
557:     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;
558:     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;
559:   }
560:   PetscLogFlops(114);
561:   VecRestoreArrayRead(Xref,  &ref);
562:   VecRestoreArray(Xreal, &real);
563:   return(0);
564: }

568: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
569: {
570:   const PetscScalar *vertices = (const PetscScalar*) ctx;
571:   const PetscScalar x0        = vertices[0];
572:   const PetscScalar y0        = vertices[1];
573:   const PetscScalar z0        = vertices[2];
574:   const PetscScalar x1        = vertices[9];
575:   const PetscScalar y1        = vertices[10];
576:   const PetscScalar z1        = vertices[11];
577:   const PetscScalar x2        = vertices[6];
578:   const PetscScalar y2        = vertices[7];
579:   const PetscScalar z2        = vertices[8];
580:   const PetscScalar x3        = vertices[3];
581:   const PetscScalar y3        = vertices[4];
582:   const PetscScalar z3        = vertices[5];
583:   const PetscScalar x4        = vertices[12];
584:   const PetscScalar y4        = vertices[13];
585:   const PetscScalar z4        = vertices[14];
586:   const PetscScalar x5        = vertices[15];
587:   const PetscScalar y5        = vertices[16];
588:   const PetscScalar z5        = vertices[17];
589:   const PetscScalar x6        = vertices[18];
590:   const PetscScalar y6        = vertices[19];
591:   const PetscScalar z6        = vertices[20];
592:   const PetscScalar x7        = vertices[21];
593:   const PetscScalar y7        = vertices[22];
594:   const PetscScalar z7        = vertices[23];
595:   const PetscScalar f_xy      = x2 - x1 - x3 + x0;
596:   const PetscScalar g_xy      = y2 - y1 - y3 + y0;
597:   const PetscScalar h_xy      = z2 - z1 - z3 + z0;
598:   const PetscScalar f_yz      = x7 - x3 - x4 + x0;
599:   const PetscScalar g_yz      = y7 - y3 - y4 + y0;
600:   const PetscScalar h_yz      = z7 - z3 - z4 + z0;
601:   const PetscScalar f_xz      = x5 - x1 - x4 + x0;
602:   const PetscScalar g_xz      = y5 - y1 - y4 + y0;
603:   const PetscScalar h_xz      = z5 - z1 - z4 + z0;
604:   const PetscScalar f_xyz     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
605:   const PetscScalar g_xyz     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
606:   const PetscScalar h_xyz     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
607:   const PetscScalar *ref;
608:   PetscErrorCode    ierr;

611:   VecGetArrayRead(Xref,  &ref);
612:   {
613:     const PetscScalar x       = ref[0];
614:     const PetscScalar y       = ref[1];
615:     const PetscScalar z       = ref[2];
616:     const PetscInt    rows[3] = {0, 1, 2};
617:     PetscScalar       values[9];

619:     values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
620:     values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
621:     values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
622:     values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
623:     values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
624:     values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
625:     values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
626:     values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
627:     values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;

629:     MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
630:   }
631:   PetscLogFlops(152);
632:   VecRestoreArrayRead(Xref,  &ref);
633:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
634:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
635:   return(0);
636: }

640: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
641: {
642:   DM             dmCoord;
643:   SNES           snes;
644:   KSP            ksp;
645:   PC             pc;
646:   Vec            coordsLocal, r, ref, real;
647:   Mat            J;
648:   const PetscScalar *coords;
649:   PetscScalar    *a;
650:   PetscInt       p;

654:   DMGetCoordinatesLocal(dm, &coordsLocal);
655:   DMGetCoordinateDM(dm, &dmCoord);
656:   SNESCreate(PETSC_COMM_SELF, &snes);
657:   SNESSetOptionsPrefix(snes, "hex_interp_");
658:   VecCreate(PETSC_COMM_SELF, &r);
659:   VecSetSizes(r, 3, 3);
660:   VecSetType(r,dm->vectype);
661:   VecDuplicate(r, &ref);
662:   VecDuplicate(r, &real);
663:   MatCreate(PETSC_COMM_SELF, &J);
664:   MatSetSizes(J, 3, 3, 3, 3);
665:   MatSetType(J, MATSEQDENSE);
666:   MatSetUp(J);
667:   SNESSetFunction(snes, r, HexMap_Private, NULL);
668:   SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
669:   SNESGetKSP(snes, &ksp);
670:   KSPGetPC(ksp, &pc);
671:   PCSetType(pc, PCLU);
672:   SNESSetFromOptions(snes);

674:   VecGetArrayRead(ctx->coords, &coords);
675:   VecGetArray(v, &a);
676:   for (p = 0; p < ctx->n; ++p) {
677:     PetscScalar *x = NULL, *vertices = NULL;
678:     PetscScalar *xi;
679:     PetscReal    xir[3];
680:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

682:     /* Can make this do all points at once */
683:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
684:     if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);
685:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
686:     if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);
687:     SNESSetFunction(snes, NULL, NULL, (void*) vertices);
688:     SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
689:     VecGetArray(real, &xi);
690:     xi[0]  = coords[p*ctx->dim+0];
691:     xi[1]  = coords[p*ctx->dim+1];
692:     xi[2]  = coords[p*ctx->dim+2];
693:     VecRestoreArray(real, &xi);
694:     SNESSolve(snes, real, ref);
695:     VecGetArray(ref, &xi);
696:     xir[0] = PetscRealPart(xi[0]);
697:     xir[1] = PetscRealPart(xi[1]);
698:     xir[2] = PetscRealPart(xi[2]);
699:     for (comp = 0; comp < ctx->dof; ++comp) {
700:       a[p*ctx->dof+comp] =
701:         x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
702:         x[3*ctx->dof+comp]*    xir[0]*(1-xir[1])*(1-xir[2]) +
703:         x[2*ctx->dof+comp]*    xir[0]*    xir[1]*(1-xir[2]) +
704:         x[1*ctx->dof+comp]*(1-xir[0])*    xir[1]*(1-xir[2]) +
705:         x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*   xir[2] +
706:         x[5*ctx->dof+comp]*    xir[0]*(1-xir[1])*   xir[2] +
707:         x[6*ctx->dof+comp]*    xir[0]*    xir[1]*   xir[2] +
708:         x[7*ctx->dof+comp]*(1-xir[0])*    xir[1]*   xir[2];
709:     }
710:     VecRestoreArray(ref, &xi);
711:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
712:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
713:   }
714:   VecRestoreArray(v, &a);
715:   VecRestoreArrayRead(ctx->coords, &coords);

717:   SNESDestroy(&snes);
718:   VecDestroy(&r);
719:   VecDestroy(&ref);
720:   VecDestroy(&real);
721:   MatDestroy(&J);
722:   return(0);
723: }

727: /*
728:   Input Parameters:
729: + ctx - The DMInterpolationInfo context
730: . dm  - The DM
731: - x   - The local vector containing the field to be interpolated

733:   Output Parameters:
734: . v   - The vector containing the interpolated values
735: */
736: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
737: {
738:   PetscInt       dim, coneSize, n;

745:   VecGetLocalSize(v, &n);
746:   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);
747:   if (n) {
748:     DMGetDimension(dm, &dim);
749:     DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);
750:     if (dim == 2) {
751:       if (coneSize == 3) {
752:         DMInterpolate_Triangle_Private(ctx, dm, x, v);
753:       } else if (coneSize == 4) {
754:         DMInterpolate_Quad_Private(ctx, dm, x, v);
755:       } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
756:     } else if (dim == 3) {
757:       if (coneSize == 4) {
758:         DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
759:       } else {
760:         DMInterpolate_Hex_Private(ctx, dm, x, v);
761:       }
762:     } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
763:   }
764:   return(0);
765: }

769: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
770: {

775:   VecDestroy(&(*ctx)->coords);
776:   PetscFree((*ctx)->points);
777:   PetscFree((*ctx)->cells);
778:   PetscFree(*ctx);
779:   *ctx = NULL;
780:   return(0);
781: }

785: /*@C
786:   SNESMonitorFields - Monitors the residual for each field separately

788:   Collective on SNES

790:   Input Parameters:
791: + snes   - the SNES context
792: . its    - iteration number
793: . fgnorm - 2-norm of residual
794: - vf  - PetscViewerAndFormat of type ASCII

796:   Notes:
797:   This routine prints the residual norm at each iteration.

799:   Level: intermediate

801: .keywords: SNES, nonlinear, default, monitor, norm
802: .seealso: SNESMonitorSet(), SNESMonitorDefault()
803: @*/
804: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
805: {
806:   PetscViewer        viewer = vf->viewer;
807:   Vec                res;
808:   DM                 dm;
809:   PetscSection       s;
810:   const PetscScalar *r;
811:   PetscReal         *lnorms, *norms;
812:   PetscInt           numFields, f, pStart, pEnd, p;
813:   PetscErrorCode     ierr;

817:   SNESGetFunction(snes, &res, 0, 0);
818:   SNESGetDM(snes, &dm);
819:   DMGetDefaultSection(dm, &s);
820:   PetscSectionGetNumFields(s, &numFields);
821:   PetscSectionGetChart(s, &pStart, &pEnd);
822:   PetscCalloc2(numFields, &lnorms, numFields, &norms);
823:   VecGetArrayRead(res, &r);
824:   for (p = pStart; p < pEnd; ++p) {
825:     for (f = 0; f < numFields; ++f) {
826:       PetscInt fdof, foff, d;

828:       PetscSectionGetFieldDof(s, p, f, &fdof);
829:       PetscSectionGetFieldOffset(s, p, f, &foff);
830:       for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
831:     }
832:   }
833:   VecRestoreArrayRead(res, &r);
834:   MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
835:   PetscViewerPushFormat(viewer,vf->format);
836:   PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
837:   PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
838:   for (f = 0; f < numFields; ++f) {
839:     if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
840:     PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
841:   }
842:   PetscViewerASCIIPrintf(viewer, "]\n");
843:   PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
844:   PetscViewerPopFormat(viewer);
845:   PetscFree2(lnorms, norms);
846:   return(0);
847: }

849: /********************* Residual Computation **************************/

853: /*@
854:   DMPlexSNESGetGeometryFEM - Return precomputed geometric data

856:   Input Parameter:
857: . dm - The DM

859:   Output Parameters:
860: . cellgeom - The values precomputed from cell geometry

862:   Level: developer

864: .seealso: DMPlexSNESSetFunctionLocal()
865: @*/
866: PetscErrorCode DMPlexSNESGetGeometryFEM(DM dm, Vec *cellgeom)
867: {
868:   DMSNES         dmsnes;
869:   PetscObject    obj;

874:   DMGetDMSNES(dm, &dmsnes);
875:   PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", &obj);
876:   if (!obj) {
877:     Vec cellgeom;

879:     DMPlexComputeGeometryFEM(dm, &cellgeom);
880:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", (PetscObject) cellgeom);
881:     VecDestroy(&cellgeom);
882:   }
884:   return(0);
885: }

889: /*@
890:   DMPlexSNESGetGeometryFVM - Return precomputed geometric data

892:   Input Parameter:
893: . dm - The DM

895:   Output Parameters:
896: + facegeom - The values precomputed from face geometry
897: . cellgeom - The values precomputed from cell geometry
898: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell

900:   Level: developer

902: .seealso: DMPlexTSSetRHSFunctionLocal()
903: @*/
904: PetscErrorCode DMPlexSNESGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
905: {
906:   DMSNES         dmsnes;
907:   PetscObject    obj;

912:   DMGetDMSNES(dm, &dmsnes);
913:   PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_facegeom_fvm", &obj);
914:   if (!obj) {
915:     Vec cellgeom, facegeom;

917:     DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);
918:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_facegeom_fvm", (PetscObject) facegeom);
919:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fvm", (PetscObject) cellgeom);
920:     VecDestroy(&facegeom);
921:     VecDestroy(&cellgeom);
922:   }
925:   if (minRadius) {DMPlexGetMinRadius(dm, minRadius);}
926:   return(0);
927: }

931: /*@
932:   DMPlexSNESGetGradientDM - Return gradient data layout

934:   Input Parameters:
935: + dm - The DM
936: - fv - The PetscFV

938:   Output Parameter:
939: . dmGrad - The layout for gradient values

941:   Level: developer

943: .seealso: DMPlexSNESGetGeometryFVM()
944: @*/
945: PetscErrorCode DMPlexSNESGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
946: {
947:   DMSNES         dmsnes;
948:   PetscObject    obj;
949:   PetscBool      computeGradients;

956:   PetscFVGetComputeGradients(fv, &computeGradients);
957:   if (!computeGradients) {*dmGrad = NULL; return(0);}
958:   DMGetDMSNES(dm, &dmsnes);
959:   PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", &obj);
960:   if (!obj) {
961:     DM  dmGrad;
962:     Vec faceGeometry, cellGeometry;

964:     DMPlexSNESGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);
965:     DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);
966:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject) dmGrad);
967:     DMDestroy(&dmGrad);
968:   }
969:   PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject *) dmGrad);
970:   return(0);
971: }

975: /*@C
976:   DMPlexGetCellFields - Retrieve the field values values for a chunk of cells

978:   Input Parameters:
979: + dm     - The DM
980: . cStart - The first cell to include
981: . cEnd   - The first cell to exclude
982: . locX   - A local vector with the solution fields
983: . locX_t - A local vector with solution field time derivatives, or NULL
984: - locA   - A local vector with auxiliary fields, or NULL

986:   Output Parameters:
987: + u   - The field coefficients
988: . u_t - The fields derivative coefficients
989: - a   - The auxiliary field coefficients

991:   Level: developer

993: .seealso: DMPlexGetFaceFields()
994: @*/
995: PetscErrorCode DMPlexGetCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
996: {
997:   DM             dmAux;
998:   PetscSection   section, sectionAux;
999:   PetscDS        prob;
1000:   PetscInt       numCells = cEnd - cStart, totDim, totDimAux, c;

1011:   DMGetDefaultSection(dm, &section);
1012:   DMGetDS(dm, &prob);
1013:   PetscDSGetTotalDimension(prob, &totDim);
1014:   if (locA) {
1015:     PetscDS probAux;

1017:     VecGetDM(locA, &dmAux);
1018:     DMGetDefaultSection(dmAux, &sectionAux);
1019:     DMGetDS(dmAux, &probAux);
1020:     PetscDSGetTotalDimension(probAux, &totDimAux);
1021:   }
1022:   DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u);
1023:   if (locX_t) {DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u_t);} else {*u_t = NULL;}
1024:   if (locA)   {DMGetWorkArray(dm, numCells*totDimAux, PETSC_SCALAR, a);} else {*a = NULL;}
1025:   for (c = cStart; c < cEnd; ++c) {
1026:     PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
1027:     PetscInt     i;

1029:     DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1030:     for (i = 0; i < totDim; ++i) ul[(c-cStart)*totDim+i] = x[i];
1031:     DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1032:     if (locX_t) {
1033:       DMPlexVecGetClosure(dm, section, locX_t, c, NULL, &x_t);
1034:       for (i = 0; i < totDim; ++i) ul_t[(c-cStart)*totDim+i] = x_t[i];
1035:       DMPlexVecRestoreClosure(dm, section, locX_t, c, NULL, &x_t);
1036:     }
1037:     if (locA) {
1038:       DM dmAuxPlex;

1040:       DMSNESConvertPlex(dmAux, &dmAuxPlex, PETSC_FALSE);
1041:       DMPlexVecGetClosure(dmAuxPlex, sectionAux, locA, c, NULL, &x);
1042:       for (i = 0; i < totDimAux; ++i) al[(c-cStart)*totDimAux+i] = x[i];
1043:       DMPlexVecRestoreClosure(dmAuxPlex, sectionAux, locA, c, NULL, &x);
1044:       DMDestroy(&dmAuxPlex);
1045:     }
1046:   }
1047:   return(0);
1048: }

1052: /*@C
1053:   DMPlexRestoreCellFields - Restore the field values values for a chunk of cells

1055:   Input Parameters:
1056: + dm     - The DM
1057: . cStart - The first cell to include
1058: . cEnd   - The first cell to exclude
1059: . locX   - A local vector with the solution fields
1060: . locX_t - A local vector with solution field time derivatives, or NULL
1061: - locA   - A local vector with auxiliary fields, or NULL

1063:   Output Parameters:
1064: + u   - The field coefficients
1065: . u_t - The fields derivative coefficients
1066: - a   - The auxiliary field coefficients

1068:   Level: developer

1070: .seealso: DMPlexGetFaceFields()
1071: @*/
1072: PetscErrorCode DMPlexRestoreCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
1073: {

1077:   DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u);
1078:   if (*u_t) {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u_t);}
1079:   if (*a)   {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, a);}
1080:   return(0);
1081: }

1085: /*@C
1086:   DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces

1088:   Input Parameters:
1089: + dm     - The DM
1090: . fStart - The first face to include
1091: . fEnd   - The first face to exclude
1092: . locX   - A local vector with the solution fields
1093: . locX_t - A local vector with solution field time derivatives, or NULL
1094: . faceGeometry - A local vector with face geometry
1095: . cellGeometry - A local vector with cell geometry
1096: - locaGrad - A local vector with field gradients, or NULL

1098:   Output Parameters:
1099: + uL - The field values at the left side of the face
1100: - uR - The field values at the right side of the face

1102:   Level: developer

1104: .seealso: DMPlexGetCellFields()
1105: @*/
1106: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1107: {
1108:   DM                 dmFace, dmCell, dmGrad = NULL;
1109:   PetscSection       section;
1110:   PetscDS            prob;
1111:   DMLabel            ghostLabel;
1112:   const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
1113:   PetscBool         *isFE;
1114:   PetscInt           dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
1115:   PetscErrorCode     ierr;

1126:   DMGetDimension(dm, &dim);
1127:   DMGetDS(dm, &prob);
1128:   DMGetDefaultSection(dm, &section);
1129:   PetscDSGetNumFields(prob, &Nf);
1130:   PetscDSGetTotalComponents(prob, &Nc);
1131:   PetscMalloc1(Nf, &isFE);
1132:   for (f = 0; f < Nf; ++f) {
1133:     PetscObject  obj;
1134:     PetscClassId id;

1136:     DMGetField(dm, f, &obj);
1137:     PetscObjectGetClassId(obj, &id);
1138:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
1139:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
1140:     else                            {isFE[f] = PETSC_FALSE;}
1141:   }
1142:   DMGetLabel(dm, "ghost", &ghostLabel);
1143:   VecGetArrayRead(locX, &x);
1144:   VecGetDM(faceGeometry, &dmFace);
1145:   VecGetArrayRead(faceGeometry, &facegeom);
1146:   VecGetDM(cellGeometry, &dmCell);
1147:   VecGetArrayRead(cellGeometry, &cellgeom);
1148:   if (locGrad) {
1149:     VecGetDM(locGrad, &dmGrad);
1150:     VecGetArrayRead(locGrad, &lgrad);
1151:   }
1152:   DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uL);
1153:   DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uR);
1154:   /* Right now just eat the extra work for FE (could make a cell loop) */
1155:   for (face = fStart, iface = 0; face < fEnd; ++face) {
1156:     const PetscInt        *cells;
1157:     PetscFVFaceGeom       *fg;
1158:     PetscFVCellGeom       *cgL, *cgR;
1159:     PetscScalar           *xL, *xR, *gL, *gR;
1160:     PetscScalar           *uLl = *uL, *uRl = *uR;
1161:     PetscInt               ghost, nsupp;

1163:     DMLabelGetValue(ghostLabel, face, &ghost);
1164:     DMPlexGetSupportSize(dm, face, &nsupp);
1165:     if (ghost >= 0) continue;
1166:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1167:     DMPlexGetSupport(dm, face, &cells);
1168:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1169:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1170:     for (f = 0; f < Nf; ++f) {
1171:       PetscInt off;

1173:       PetscDSGetComponentOffset(prob, f, &off);
1174:       if (isFE[f]) {
1175:         const PetscInt *cone;
1176:         PetscInt        comp, coneSize, faceLocL, faceLocR, ldof, rdof, d;

1178:         xL = xR = NULL;
1179:         DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1180:         DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1181:         DMPlexGetCone(dm, cells[0], &cone);
1182:         DMPlexGetConeSize(dm, cells[0], &coneSize);
1183:         for (faceLocL = 0; faceLocL < coneSize; ++faceLocL) if (cone[faceLocL] == face) break;
1184:         if (faceLocL == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[0]);
1185:         DMPlexGetCone(dm, cells[1], &cone);
1186:         DMPlexGetConeSize(dm, cells[1], &coneSize);
1187:         for (faceLocR = 0; faceLocR < coneSize; ++faceLocR) if (cone[faceLocR] == face) break;
1188:         if (faceLocR == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[1]);
1189:         /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
1190:         EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
1191:         if (rdof == ldof) {EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
1192:         else              {PetscSectionGetFieldComponents(section, f, &comp); for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
1193:         DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1194:         DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1195:       } else {
1196:         PetscFV  fv;
1197:         PetscInt numComp, c;

1199:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
1200:         PetscFVGetNumComponents(fv, &numComp);
1201:         if (nsupp > 2) {
1202:           for (f = 0; f < Nf; ++f) {
1203:             PetscInt off;

1205:             PetscDSGetComponentOffset(prob, f, &off);
1206:             PetscFVGetNumComponents(fv, &numComp);
1207:             for (c = 0; c < numComp; ++c) {
1208:               uLl[iface*Nc+off+c] = 0.;
1209:               uRl[iface*Nc+off+c] = 0.;
1210:             }
1211:           }
1212:           continue;
1213:         }
1214:         DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
1215:         DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
1216:         if (dmGrad) {
1217:           PetscReal dxL[3], dxR[3];

1219:           DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
1220:           DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
1221:           DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
1222:           DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
1223:           for (c = 0; c < numComp; ++c) {
1224:             uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
1225:             uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
1226:           }
1227:         } else {
1228:           for (c = 0; c < numComp; ++c) {
1229:             uLl[iface*Nc+off+c] = xL[c];
1230:             uRl[iface*Nc+off+c] = xR[c];
1231:           }
1232:         }
1233:       }
1234:     }
1235:     ++iface;
1236:   }
1237:   VecRestoreArrayRead(locX, &x);
1238:   VecRestoreArrayRead(faceGeometry, &facegeom);
1239:   VecRestoreArrayRead(cellGeometry, &cellgeom);
1240:   if (locGrad) {
1241:     VecRestoreArrayRead(locGrad, &lgrad);
1242:   }
1243:   PetscFree(isFE);
1244:   return(0);
1245: }

1249: /*@C
1250:   DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces

1252:   Input Parameters:
1253: + dm     - The DM
1254: . fStart - The first face to include
1255: . fEnd   - The first face to exclude
1256: . locX   - A local vector with the solution fields
1257: . locX_t - A local vector with solution field time derivatives, or NULL
1258: . faceGeometry - A local vector with face geometry
1259: . cellGeometry - A local vector with cell geometry
1260: - locaGrad - A local vector with field gradients, or NULL

1262:   Output Parameters:
1263: + uL - The field values at the left side of the face
1264: - uR - The field values at the right side of the face

1266:   Level: developer

1268: .seealso: DMPlexGetFaceFields()
1269: @*/
1270: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1271: {

1275:   DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uL);
1276:   DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uR);
1277:   return(0);
1278: }

1282: /*@C
1283:   DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces

1285:   Input Parameters:
1286: + dm     - The DM
1287: . fStart - The first face to include
1288: . fEnd   - The first face to exclude
1289: . faceGeometry - A local vector with face geometry
1290: - cellGeometry - A local vector with cell geometry

1292:   Output Parameters:
1293: + fgeom - The extract the face centroid and normal
1294: - vol   - The cell volume

1296:   Level: developer

1298: .seealso: DMPlexGetCellFields()
1299: @*/
1300: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1301: {
1302:   DM                 dmFace, dmCell;
1303:   DMLabel            ghostLabel;
1304:   const PetscScalar *facegeom, *cellgeom;
1305:   PetscInt           dim, numFaces = fEnd - fStart, iface, face;
1306:   PetscErrorCode     ierr;

1314:   DMGetDimension(dm, &dim);
1315:   DMGetLabel(dm, "ghost", &ghostLabel);
1316:   VecGetDM(faceGeometry, &dmFace);
1317:   VecGetArrayRead(faceGeometry, &facegeom);
1318:   VecGetDM(cellGeometry, &dmCell);
1319:   VecGetArrayRead(cellGeometry, &cellgeom);
1320:   PetscMalloc1(numFaces, fgeom);
1321:   DMGetWorkArray(dm, numFaces*2, PETSC_SCALAR, vol);
1322:   for (face = fStart, iface = 0; face < fEnd; ++face) {
1323:     const PetscInt        *cells;
1324:     PetscFVFaceGeom       *fg;
1325:     PetscFVCellGeom       *cgL, *cgR;
1326:     PetscFVFaceGeom       *fgeoml = *fgeom;
1327:     PetscReal             *voll   = *vol;
1328:     PetscInt               ghost, d;

1330:     DMLabelGetValue(ghostLabel, face, &ghost);
1331:     if (ghost >= 0) continue;
1332:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1333:     DMPlexGetSupport(dm, face, &cells);
1334:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1335:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1336:     for (d = 0; d < dim; ++d) {
1337:       fgeoml[iface].centroid[d] = fg->centroid[d];
1338:       fgeoml[iface].normal[d]   = fg->normal[d];
1339:     }
1340:     voll[iface*2+0] = cgL->volume;
1341:     voll[iface*2+1] = cgR->volume;
1342:     ++iface;
1343:   }
1344:   VecRestoreArrayRead(faceGeometry, &facegeom);
1345:   VecRestoreArrayRead(cellGeometry, &cellgeom);
1346:   return(0);
1347: }

1351: /*@C
1352:   DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces

1354:   Input Parameters:
1355: + dm     - The DM
1356: . fStart - The first face to include
1357: . fEnd   - The first face to exclude
1358: . faceGeometry - A local vector with face geometry
1359: - cellGeometry - A local vector with cell geometry

1361:   Output Parameters:
1362: + fgeom - The extract the face centroid and normal
1363: - vol   - The cell volume

1365:   Level: developer

1367: .seealso: DMPlexGetFaceFields()
1368: @*/
1369: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1370: {

1374:   PetscFree(*fgeom);
1375:   DMRestoreWorkArray(dm, 0, PETSC_REAL, vol);
1376:   return(0);
1377: }

1381: static PetscErrorCode DMPlexApplyLimiter_Internal (DM dm, DM dmCell, PetscLimiter lim, PetscInt dim, PetscInt totDim, PetscInt cell, PetscInt face, PetscInt fStart, PetscInt fEnd, PetscReal *cellPhi, const PetscScalar *x,
1382:                                                    const PetscScalar *cellgeom, const PetscFVCellGeom *cg, const PetscScalar *cx, const PetscScalar *cgrad)
1383: {
1384:   const PetscInt        *children;
1385:   PetscInt               numChildren;
1386:   PetscErrorCode         ierr;

1389:   DMPlexGetTreeChildren(dm,face,&numChildren,&children);
1390:   if (numChildren) {
1391:     PetscInt c;

1393:     for (c = 0; c < numChildren; c++) {
1394:       PetscInt childFace = children[c];

1396:       if (childFace >= fStart && childFace < fEnd) {
1397:         DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,totDim,cell,childFace,fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
1398:       }
1399:     }
1400:   }
1401:   else {
1402:     PetscScalar           *ncx;
1403:     PetscFVCellGeom       *ncg;
1404:     const PetscInt        *fcells;
1405:     PetscInt               ncell, d;
1406:     PetscReal              v[3];

1408:     DMPlexGetSupport(dm, face, &fcells);
1409:     ncell = cell == fcells[0] ? fcells[1] : fcells[0];
1410:     DMPlexPointLocalRead(dm, ncell, x, &ncx);
1411:     DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);
1412:     DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v);
1413:     for (d = 0; d < totDim; ++d) {
1414:       /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
1415:       PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / DMPlex_DotD_Internal(dim, &cgrad[d*dim], v);

1417:       PetscLimiterLimit(lim, flim, &phi);
1418:       cellPhi[d] = PetscMin(cellPhi[d], phi);
1419:     }
1420:   }
1421:   return(0);
1422: }

1426: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
1427: {
1428:   DM                 dmFace, dmCell, dmGrad;
1429:   DMLabel            ghostLabel;
1430:   PetscDS            prob;
1431:   PetscFV            fvm;
1432:   PetscLimiter       lim;
1433:   const PetscScalar *facegeom, *cellgeom, *x;
1434:   PetscScalar       *gr;
1435:   PetscReal         *cellPhi;
1436:   PetscInt           dim, face, cell, totDim, cStart, cEnd, cEndInterior;
1437:   PetscErrorCode     ierr;

1440:   DMGetDimension(dm, &dim);
1441:   DMGetDS(dm, &prob);
1442:   PetscDSGetTotalDimension(prob, &totDim);
1443:   DMGetLabel(dm, "ghost", &ghostLabel);
1444:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fvm);
1445:   PetscFVGetLimiter(fvm, &lim);
1446:   VecGetDM(faceGeometry, &dmFace);
1447:   VecGetArrayRead(faceGeometry, &facegeom);
1448:   VecGetDM(cellGeometry, &dmCell);
1449:   VecGetArrayRead(cellGeometry, &cellgeom);
1450:   VecGetArrayRead(locX, &x);
1451:   VecGetDM(grad, &dmGrad);
1452:   VecZeroEntries(grad);
1453:   VecGetArray(grad, &gr);
1454:   /* Reconstruct gradients */
1455:   for (face = fStart; face < fEnd; ++face) {
1456:     const PetscInt        *cells;
1457:     PetscFVFaceGeom       *fg;
1458:     PetscScalar           *cx[2];
1459:     PetscScalar           *cgrad[2];
1460:     PetscBool              boundary;
1461:     PetscInt               ghost, c, pd, d, numChildren, numCells;

1463:     DMLabelGetValue(ghostLabel, face, &ghost);
1464:     DMIsBoundaryPoint(dm, face, &boundary);
1465:     DMPlexGetTreeChildren(dm, face, &numChildren, NULL);
1466:     if (ghost >= 0 || boundary || numChildren) continue;
1467:     DMPlexGetSupportSize(dm, face, &numCells);
1468:     if (numCells != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "facet %d has %d support points: expected 2",face,numCells);
1469:     DMPlexGetSupport(dm, face, &cells);
1470:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1471:     for (c = 0; c < 2; ++c) {
1472:       DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
1473:       DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]);
1474:     }
1475:     for (pd = 0; pd < totDim; ++pd) {
1476:       PetscScalar delta = cx[1][pd] - cx[0][pd];

1478:       for (d = 0; d < dim; ++d) {
1479:         if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta;
1480:         if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta;
1481:       }
1482:     }
1483:   }
1484:   /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
1485:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1486:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1487:   cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
1488:   DMGetWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1489:   for (cell = dmGrad && lim ? cStart : cEnd; cell < cEndInterior; ++cell) {
1490:     const PetscInt        *faces;
1491:     PetscScalar           *cx;
1492:     PetscFVCellGeom       *cg;
1493:     PetscScalar           *cgrad;
1494:     PetscInt               coneSize, f, pd, d;

1496:     DMPlexGetConeSize(dm, cell, &coneSize);
1497:     DMPlexGetCone(dm, cell, &faces);
1498:     DMPlexPointLocalRead(dm, cell, x, &cx);
1499:     DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
1500:     DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
1501:     if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
1502:     /* Limiter will be minimum value over all neighbors */
1503:     for (d = 0; d < totDim; ++d) cellPhi[d] = PETSC_MAX_REAL;
1504:     for (f = 0; f < coneSize; ++f) {
1505:       DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,totDim,cell,faces[f],fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
1506:     }
1507:     /* Apply limiter to gradient */
1508:     for (pd = 0; pd < totDim; ++pd)
1509:       /* Scalar limiter applied to each component separately */
1510:       for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd];
1511:   }
1512:   DMRestoreWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1513:   VecRestoreArrayRead(faceGeometry, &facegeom);
1514:   VecRestoreArrayRead(cellGeometry, &cellgeom);
1515:   VecRestoreArrayRead(locX, &x);
1516:   VecRestoreArray(grad, &gr);
1517:   return(0);
1518: }

1522: PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX_t, Vec locF, void *user)
1523: {
1524:   DM_Plex         *mesh = (DM_Plex *) dm->data;
1525:   PetscSection     section;
1526:   PetscDS          prob;
1527:   DMLabel          depth;
1528:   PetscFECellGeom *cgeom;
1529:   PetscScalar     *u = NULL, *u_t = NULL, *elemVec = NULL;
1530:   PetscInt         dim, Nf, f, totDimBd, numBd, bd;
1531:   PetscErrorCode   ierr;

1534:   DMGetDimension(dm, &dim);
1535:   DMGetDefaultSection(dm, &section);
1536:   DMGetDS(dm, &prob);
1537:   PetscDSGetNumFields(prob, &Nf);
1538:   PetscDSGetTotalBdDimension(prob, &totDimBd);
1539:   DMPlexGetDepthLabel(dm, &depth);
1540:   DMGetNumBoundary(dm, &numBd);
1541:   for (bd = 0; bd < numBd; ++bd) {
1542:     const char     *bdLabel;
1543:     DMLabel         label;
1544:     IS              pointIS;
1545:     const PetscInt *points;
1546:     const PetscInt *values;
1547:     PetscInt        field, numValues, v, numPoints, p, dep, numFaces;
1548:     PetscBool       isEssential;
1549:     PetscObject     obj;
1550:     PetscClassId    id;

1552:     DMGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
1553:     DMGetField(dm, field, &obj);
1554:     PetscObjectGetClassId(obj, &id);
1555:     if ((id != PETSCFE_CLASSID) || isEssential) continue;
1556:     DMGetLabel(dm, bdLabel, &label);
1557:     for (v = 0; v < numValues; ++v) {
1558:       DMLabelGetStratumSize(label, values[v], &numPoints);
1559:       DMLabelGetStratumIS(label, values[v], &pointIS);
1560:       if (!pointIS) continue; /* No points with that id on this process */
1561:       ISGetIndices(pointIS, &points);
1562:       for (p = 0, numFaces = 0; p < numPoints; ++p) {
1563:         DMLabelGetValue(depth, points[p], &dep);
1564:         if (dep == dim-1) ++numFaces;
1565:       }
1566:       PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd,&elemVec);
1567:       if (locX_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
1568:       for (p = 0, f = 0; p < numPoints; ++p) {
1569:         const PetscInt point = points[p];
1570:         PetscScalar   *x     = NULL;
1571:         PetscInt       i;

1573:         DMLabelGetValue(depth, points[p], &dep);
1574:         if (dep != dim-1) continue;
1575:         DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
1576:         DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
1577:         if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
1578:         /* TODO: Matt, this is wrong if feBd does not match fe: i.e., if the order differs. */
1579:         DMPlexVecGetClosure(dm, section, locX, point, NULL, &x);
1580:         for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1581:         DMPlexVecRestoreClosure(dm, section, locX, point, NULL, &x);
1582:         if (locX_t) {
1583:           DMPlexVecGetClosure(dm, section, locX_t, point, NULL, &x);
1584:           for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1585:           DMPlexVecRestoreClosure(dm, section, locX_t, point, NULL, &x);
1586:         }
1587:         ++f;
1588:       }
1589:       for (f = 0; f < Nf; ++f) {
1590:         PetscFE         fe;
1591:         PetscQuadrature q;
1592:         PetscInt        numQuadPoints, Nb;
1593:         /* Conforming batches */
1594:         PetscInt        numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1595:         /* Remainder */
1596:         PetscInt        Nr, offset;

1598:         PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);
1599:         PetscFEGetQuadrature(fe, &q);
1600:         PetscFEGetDimension(fe, &Nb);
1601:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1602:         PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1603:         blockSize = Nb*numQuadPoints;
1604:         batchSize = numBlocks * blockSize;
1605:          PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1606:         numChunks = numFaces / (numBatches*batchSize);
1607:         Ne        = numChunks*numBatches*batchSize;
1608:         Nr        = numFaces % (numBatches*batchSize);
1609:         offset    = numFaces - Nr;
1610:         PetscFEIntegrateBdResidual(fe, prob, f, Ne, cgeom, u, u_t, NULL, NULL, elemVec);
1611:         PetscFEIntegrateBdResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);
1612:       }
1613:       for (p = 0, f = 0; p < numPoints; ++p) {
1614:         const PetscInt point = points[p];

1616:         DMLabelGetValue(depth, point, &dep);
1617:         if (dep != dim-1) continue;
1618:         if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);}
1619:         DMPlexVecSetClosure(dm, NULL, locF, point, &elemVec[f*totDimBd], ADD_ALL_VALUES);
1620:         ++f;
1621:       }
1622:       ISRestoreIndices(pointIS, &points);
1623:       ISDestroy(&pointIS);
1624:       PetscFree3(u,cgeom,elemVec);
1625:       if (locX_t) {PetscFree(u_t);}
1626:     }
1627:   }
1628:   return(0);
1629: }

1633: /*@
1634:   DMPlexReconstructGradientsFVM - reconstruct the gradient of a vector using a finite volume method.

1636:   Input Parameters:
1637: + dm - the mesh
1638: - locX - the local representation of the vector

1640:   Output Parameter:
1641: . grad - the global representation of the gradient

1643:   Level: developer

1645: .seealso: DMPlexSNESGetGradientDM()
1646: @*/
1647: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
1648: {
1649:   PetscDS          prob;
1650:   PetscInt         Nf, f, fStart, fEnd;
1651:   PetscBool        useFVM = PETSC_FALSE;
1652:   PetscFV          fvm = NULL;
1653:   Vec              faceGeometryFVM, cellGeometryFVM;
1654:   PetscFVCellGeom  *cgeomFVM   = NULL;
1655:   PetscFVFaceGeom  *fgeomFVM   = NULL;
1656:   DM               dmGrad = NULL;
1657:   PetscErrorCode   ierr;

1660:   DMGetDS(dm, &prob);
1661:   PetscDSGetNumFields(prob, &Nf);
1662:   for (f = 0; f < Nf; ++f) {
1663:     PetscObject  obj;
1664:     PetscClassId id;

1666:     PetscDSGetDiscretization(prob, f, &obj);
1667:     PetscObjectGetClassId(obj, &id);
1668:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1669:   }
1670:   if (!useFVM) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm does not have a finite volume discretization");
1671:   DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1672:   if (!dmGrad) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm's finite volume discretization does not reconstruct gradients");
1673:   DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1674:   VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1675:   VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1676:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1677:   DMPlexReconstructGradients_Internal(dm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1678:   return(0);
1679: }

1683: PetscErrorCode DMPlexComputeResidual_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
1684: {
1685:   DM_Plex          *mesh       = (DM_Plex *) dm->data;
1686:   const char       *name       = "Residual";
1687:   DM                dmAux      = NULL;
1688:   DM                dmGrad     = NULL;
1689:   DMLabel           ghostLabel = NULL;
1690:   PetscDS           prob       = NULL;
1691:   PetscDS           probAux    = NULL;
1692:   PetscSection      section    = NULL;
1693:   PetscBool         useFEM     = PETSC_FALSE;
1694:   PetscBool         useFVM     = PETSC_FALSE;
1695:   PetscBool         isImplicit = (locX_t || time == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
1696:   PetscFV           fvm        = NULL;
1697:   PetscFECellGeom  *cgeomFEM   = NULL;
1698:   PetscScalar      *cgeomScal;
1699:   PetscFVCellGeom  *cgeomFVM   = NULL;
1700:   PetscFVFaceGeom  *fgeomFVM   = NULL;
1701:   Vec               locA, cellGeometryFEM = NULL, cellGeometryFVM = NULL, faceGeometryFVM = NULL, grad, locGrad = NULL;
1702:   PetscScalar      *u, *u_t, *a, *uL, *uR;
1703:   PetscInt          Nf, f, totDim, totDimAux, numChunks, cellChunkSize, faceChunkSize, chunk, fStart, fEnd;
1704:   PetscErrorCode    ierr;

1707:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
1708:   /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
1709:   /* TODO The FVM geometry is over-manipulated. Make the precalc functions return exactly what we need */
1710:   /* FEM+FVM */
1711:   /* 1: Get sizes from dm and dmAux */
1712:   DMGetDefaultSection(dm, &section);
1713:   DMGetLabel(dm, "ghost", &ghostLabel);
1714:   DMGetDS(dm, &prob);
1715:   PetscDSGetNumFields(prob, &Nf);
1716:   PetscDSGetTotalDimension(prob, &totDim);
1717:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1718:   if (locA) {
1719:     VecGetDM(locA, &dmAux);
1720:     DMGetDS(dmAux, &probAux);
1721:     PetscDSGetTotalDimension(probAux, &totDimAux);
1722:   }
1723:   /* 2: Get geometric data */
1724:   for (f = 0; f < Nf; ++f) {
1725:     PetscObject  obj;
1726:     PetscClassId id;
1727:     PetscBool    fimp;

1729:     PetscDSGetImplicit(prob, f, &fimp);
1730:     if (isImplicit != fimp) continue;
1731:     PetscDSGetDiscretization(prob, f, &obj);
1732:     PetscObjectGetClassId(obj, &id);
1733:     if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
1734:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1735:   }
1736:   if (useFEM) {
1737:     DMPlexSNESGetGeometryFEM(dm, &cellGeometryFEM);
1738:     VecGetArray(cellGeometryFEM, &cgeomScal);
1739:     if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1740:       DM dmCell;
1741:       PetscInt c;

1743:       VecGetDM(cellGeometryFEM,&dmCell);
1744:       PetscMalloc1(cEnd-cStart,&cgeomFEM);
1745:       for (c = 0; c < cEnd - cStart; c++) {
1746:         PetscScalar *thisgeom;

1748:         DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
1749:         cgeomFEM[c] = *((PetscFECellGeom *) thisgeom);
1750:       }
1751:     }
1752:     else {
1753:       cgeomFEM = (PetscFECellGeom *) cgeomScal;
1754:     }
1755:   }
1756:   if (useFVM) {
1757:     DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1758:     VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1759:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1760:     /* Reconstruct and limit cell gradients */
1761:     DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1762:     if (dmGrad) {
1763:       DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1764:       DMGetGlobalVector(dmGrad, &grad);
1765:       DMPlexReconstructGradients_Internal(dm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1766:       /* Communicate gradient values */
1767:       DMGetLocalVector(dmGrad, &locGrad);
1768:       DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1769:       DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1770:       DMRestoreGlobalVector(dmGrad, &grad);
1771:     }
1772:     /* Handle non-essential (e.g. outflow) boundary values */
1773:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, time, faceGeometryFVM, cellGeometryFVM, locGrad);
1774:   }
1775:   /* Loop over chunks */
1776:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1777:   numChunks     = 1;
1778:   cellChunkSize = (cEnd - cStart)/numChunks;
1779:   faceChunkSize = (fEnd - fStart)/numChunks;
1780:   for (chunk = 0; chunk < numChunks; ++chunk) {
1781:     PetscScalar     *elemVec, *fluxL, *fluxR;
1782:     PetscReal       *vol;
1783:     PetscFVFaceGeom *fgeom;
1784:     PetscInt         cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, cell;
1785:     PetscInt         fS = fStart+chunk*faceChunkSize, fE = PetscMin(fS+faceChunkSize, fEnd), numFaces = fE - fS, face;

1787:     /* Extract field coefficients */
1788:     if (useFEM) {
1789:       DMPlexGetCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1790:       DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1791:       PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
1792:     }
1793:     if (useFVM) {
1794:       DMPlexGetFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1795:       DMPlexGetFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1796:       DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1797:       DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1798:       PetscMemzero(fluxL, numFaces*totDim * sizeof(PetscScalar));
1799:       PetscMemzero(fluxR, numFaces*totDim * sizeof(PetscScalar));
1800:     }
1801:     /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
1802:     /* Loop over fields */
1803:     for (f = 0; f < Nf; ++f) {
1804:       PetscObject  obj;
1805:       PetscClassId id;
1806:       PetscBool    fimp;
1807:       PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1809:       PetscDSGetImplicit(prob, f, &fimp);
1810:       if (isImplicit != fimp) continue;
1811:       PetscDSGetDiscretization(prob, f, &obj);
1812:       PetscObjectGetClassId(obj, &id);
1813:       if (id == PETSCFE_CLASSID) {
1814:         PetscFE         fe = (PetscFE) obj;
1815:         PetscQuadrature q;
1816:         PetscInt        Nq, Nb;

1818:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);

1820:         PetscFEGetQuadrature(fe, &q);
1821:         PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1822:         PetscFEGetDimension(fe, &Nb);
1823:         blockSize = Nb*Nq;
1824:         batchSize = numBlocks * blockSize;
1825:          PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1826:         numChunks = numCells / (numBatches*batchSize);
1827:         Ne        = numChunks*numBatches*batchSize;
1828:         Nr        = numCells % (numBatches*batchSize);
1829:         offset    = numCells - Nr;
1830:         /* Integrate FE residual to get elemVec (need fields at quadrature points) */
1831:         /*   For FV, I think we use a P0 basis and the cell coefficients (for subdivided cells, we can tweak the basis tabulation to be the indicator function) */
1832:         PetscFEIntegrateResidual(fe, prob, f, Ne, cgeomFEM, u, u_t, probAux, a, elemVec);
1833:         PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1834:       } else if (id == PETSCFV_CLASSID) {
1835:         PetscFV fv = (PetscFV) obj;

1837:         Ne = numFaces;
1838:         /* Riemann solve over faces (need fields at face centroids) */
1839:         /*   We need to evaluate FE fields at those coordinates */
1840:         PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
1841:       } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1842:     }
1843:     if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1844:       PetscFree(cgeomFEM);
1845:     }
1846:     else {
1847:       cgeomFEM = NULL;
1848:     }
1849:     if (cellGeometryFEM) {VecRestoreArray(cellGeometryFEM, &cgeomScal);}
1850:     /* Loop over domain */
1851:     if (useFEM) {
1852:       /* Add elemVec to locX */
1853:       for (cell = cS; cell < cE; ++cell) {
1854:         if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cell*totDim]);}
1855:         DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cell*totDim], ADD_ALL_VALUES);
1856:       }
1857:     }
1858:     if (useFVM) {
1859:       PetscScalar *fa;
1860:       PetscInt     iface;

1862:       VecGetArray(locF, &fa);
1863:       for (f = 0; f < Nf; ++f) {
1864:         PetscFV      fv;
1865:         PetscObject  obj;
1866:         PetscClassId id;
1867:         PetscInt     foff, pdim;

1869:         PetscDSGetDiscretization(prob, f, &obj);
1870:         PetscDSGetFieldOffset(prob, f, &foff);
1871:         PetscObjectGetClassId(obj, &id);
1872:         if (id != PETSCFV_CLASSID) continue;
1873:         fv   = (PetscFV) obj;
1874:         PetscFVGetNumComponents(fv, &pdim);
1875:         /* Accumulate fluxes to cells */
1876:         for (face = fS, iface = 0; face < fE; ++face) {
1877:           const PetscInt *cells;
1878:           PetscScalar    *fL, *fR;
1879:           PetscInt        ghost, d, nsupp;

1881:           DMLabelGetValue(ghostLabel, face, &ghost);
1882:           DMPlexGetSupportSize(dm, face, &nsupp);
1883:           if (ghost >= 0) continue;
1884:           if (nsupp > 2) { /* noop */
1885:             ++iface;
1886:             continue;
1887:           }
1888:           DMPlexGetSupport(dm, face, &cells);
1889:           DMPlexPointGlobalFieldRef(dm, cells[0], f, fa, &fL);
1890:           DMPlexPointGlobalFieldRef(dm, cells[1], f, fa, &fR);
1891:           for (d = 0; d < pdim; ++d) {
1892:             if (fL) fL[d] -= fluxL[iface*totDim+foff+d];
1893:             if (fR) fR[d] += fluxR[iface*totDim+foff+d];
1894:           }
1895:           ++iface;
1896:         }
1897:       }
1898:       VecRestoreArray(locF, &fa);
1899:     }
1900:     /* Handle time derivative */
1901:     if (locX_t) {
1902:       PetscScalar *x_t, *fa;

1904:       VecGetArray(locF, &fa);
1905:       VecGetArray(locX_t, &x_t);
1906:       for (f = 0; f < Nf; ++f) {
1907:         PetscFV      fv;
1908:         PetscObject  obj;
1909:         PetscClassId id;
1910:         PetscInt     pdim, d;

1912:         PetscDSGetDiscretization(prob, f, &obj);
1913:         PetscObjectGetClassId(obj, &id);
1914:         if (id != PETSCFV_CLASSID) continue;
1915:         fv   = (PetscFV) obj;
1916:         PetscFVGetNumComponents(fv, &pdim);
1917:         for (cell = cS; cell < cE; ++cell) {
1918:           PetscScalar *u_t, *r;

1920:           DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
1921:           DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
1922:           for (d = 0; d < pdim; ++d) r[d] += u_t[d];
1923:         }
1924:       }
1925:       VecRestoreArray(locX_t, &x_t);
1926:       VecRestoreArray(locF, &fa);
1927:     }
1928:     if (useFEM) {
1929:       DMPlexRestoreCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1930:       DMRestoreWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1931:     }
1932:     if (useFVM) {
1933:       DMPlexRestoreFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1934:       DMPlexRestoreFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1935:       DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1936:       DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1937:       if (dmGrad) {DMRestoreLocalVector(dmGrad, &locGrad);}
1938:     }
1939:   }

1941:   if (useFEM) {DMPlexComputeBdResidual_Internal(dm, locX, locX_t, locF, user);}

1943:   /* FEM */
1944:   /* 1: Get sizes from dm and dmAux */
1945:   /* 2: Get geometric data */
1946:   /* 3: Handle boundary values */
1947:   /* 4: Loop over domain */
1948:   /*   Extract coefficients */
1949:   /* Loop over fields */
1950:   /*   Set tiling for FE*/
1951:   /*   Integrate FE residual to get elemVec */
1952:   /*     Loop over subdomain */
1953:   /*       Loop over quad points */
1954:   /*         Transform coords to real space */
1955:   /*         Evaluate field and aux fields at point */
1956:   /*         Evaluate residual at point */
1957:   /*         Transform residual to real space */
1958:   /*       Add residual to elemVec */
1959:   /* Loop over domain */
1960:   /*   Add elemVec to locX */

1962:   /* FVM */
1963:   /* Get geometric data */
1964:   /* If using gradients */
1965:   /*   Compute gradient data */
1966:   /*   Loop over domain faces */
1967:   /*     Count computational faces */
1968:   /*     Reconstruct cell gradient */
1969:   /*   Loop over domain cells */
1970:   /*     Limit cell gradients */
1971:   /* Handle boundary values */
1972:   /* Loop over domain faces */
1973:   /*   Read out field, centroid, normal, volume for each side of face */
1974:   /* Riemann solve over faces */
1975:   /* Loop over domain faces */
1976:   /*   Accumulate fluxes to cells */
1977:   /* TODO Change printFEM to printDisc here */
1978:   if (mesh->printFEM) {
1979:     Vec         locFbc;
1980:     PetscInt    pStart, pEnd, p, maxDof;
1981:     PetscScalar *zeroes;

1983:     VecDuplicate(locF,&locFbc);
1984:     VecCopy(locF,locFbc);
1985:     PetscSectionGetChart(section,&pStart,&pEnd);
1986:     PetscSectionGetMaxDof(section,&maxDof);
1987:     PetscCalloc1(maxDof,&zeroes);
1988:     for (p = pStart; p < pEnd; p++) {
1989:       VecSetValuesSection(locFbc,section,p,zeroes,INSERT_BC_VALUES);
1990:     }
1991:     PetscFree(zeroes);
1992:     DMPrintLocalVec(dm, name, mesh->printTol, locFbc);
1993:     VecDestroy(&locFbc);
1994:   }
1995:   PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
1996:   return(0);
1997: }

2001: static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
2002: {
2003:   DM                dmCh, dmAux;
2004:   Vec               A, cellgeom;
2005:   PetscDS           prob, probCh, probAux = NULL;
2006:   PetscQuadrature   q;
2007:   PetscSection      section, sectionAux;
2008:   PetscFECellGeom  *cgeom = NULL;
2009:   PetscScalar      *cgeomScal;
2010:   PetscScalar      *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
2011:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
2012:   PetscInt          totDim, totDimAux, diffCell = 0;
2013:   PetscErrorCode    ierr;

2016:   DMGetDimension(dm, &dim);
2017:   DMGetDefaultSection(dm, &section);
2018:   DMGetDS(dm, &prob);
2019:   PetscDSGetTotalDimension(prob, &totDim);
2020:   PetscSectionGetNumFields(section, &Nf);
2021:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2022:   numCells = cEnd - cStart;
2023:   PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);
2024:   DMGetDS(dmCh, &probCh);
2025:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2026:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2027:   if (dmAux) {
2028:     DMGetDefaultSection(dmAux, &sectionAux);
2029:     DMGetDS(dmAux, &probAux);
2030:     PetscDSGetTotalDimension(probAux, &totDimAux);
2031:   }
2032:   VecSet(F, 0.0);
2033:   PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);
2034:   PetscMalloc1(numCells*totDim,&elemVecCh);
2035:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2036:   DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2037:   VecGetArray(cellgeom, &cgeomScal);
2038:   if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2039:     DM dmCell;

2041:     VecGetDM(cellgeom,&dmCell);
2042:     PetscMalloc1(cEnd-cStart,&cgeom);
2043:     for (c = 0; c < cEnd - cStart; c++) {
2044:       PetscScalar *thisgeom;

2046:       DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2047:       cgeom[c] = *((PetscFECellGeom *) thisgeom);
2048:     }
2049:   }
2050:   else {
2051:     cgeom = (PetscFECellGeom *) cgeomScal;
2052:   }
2053:   for (c = cStart; c < cEnd; ++c) {
2054:     PetscScalar *x = NULL, *x_t = NULL;
2055:     PetscInt     i;

2057:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2058:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
2059:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2060:     if (X_t) {
2061:       DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2062:       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
2063:       DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2064:     }
2065:     if (dmAux) {
2066:       DM dmAuxPlex;

2068:       DMSNESConvertPlex(dmAux,&dmAuxPlex, PETSC_FALSE);
2069:       DMPlexVecGetClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
2070:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
2071:       DMPlexVecRestoreClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
2072:       DMDestroy(&dmAuxPlex);
2073:     }
2074:   }
2075:   for (f = 0; f < Nf; ++f) {
2076:     PetscFE  fe, feCh;
2077:     PetscInt numQuadPoints, Nb;
2078:     /* Conforming batches */
2079:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2080:     /* Remainder */
2081:     PetscInt Nr, offset;

2083:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
2084:     PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);
2085:     PetscFEGetQuadrature(fe, &q);
2086:     PetscFEGetDimension(fe, &Nb);
2087:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2088:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
2089:     blockSize = Nb*numQuadPoints;
2090:     batchSize = numBlocks * blockSize;
2091:      PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2092:     numChunks = numCells / (numBatches*batchSize);
2093:     Ne        = numChunks*numBatches*batchSize;
2094:     Nr        = numCells % (numBatches*batchSize);
2095:     offset    = numCells - Nr;
2096:     PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVec);
2097:     PetscFEIntegrateResidual(feCh, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVecCh);
2098:     PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
2099:     PetscFEIntegrateResidual(feCh, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVecCh[offset*totDim]);
2100:   }
2101:   for (c = cStart; c < cEnd; ++c) {
2102:     PetscBool diff = PETSC_FALSE;
2103:     PetscInt  d;

2105:     for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
2106:     if (diff) {
2107:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);
2108:       DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);
2109:       DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);
2110:       ++diffCell;
2111:     }
2112:     if (diffCell > 9) break;
2113:     DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_ALL_VALUES);
2114:   }
2115:   if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2116:     PetscFree(cgeom);
2117:   }
2118:   else {
2119:     cgeom = NULL;
2120:   }
2121:   VecRestoreArray(cellgeom, &cgeomScal);
2122:   PetscFree3(u,u_t,elemVec);
2123:   PetscFree(elemVecCh);
2124:   if (dmAux) {PetscFree(a);}
2125:   return(0);
2126: }

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

2133:   Input Parameters:
2134: + dm - The mesh
2135: . X  - Local solution
2136: - user - The user context

2138:   Output Parameter:
2139: . F  - Local output vector

2141:   Level: developer

2143: .seealso: DMPlexComputeJacobianActionFEM()
2144: @*/
2145: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
2146: {
2147:   PetscObject    check;
2148:   PetscInt       cStart, cEnd, cEndInterior;
2149:   DM             plex;

2153:   DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2154:   DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2155:   DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2156:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2157:   /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
2158:   PetscObjectQuery((PetscObject) plex, "dmCh", &check);
2159:   if (check) {DMPlexComputeResidualFEM_Check_Internal(plex, X, NULL, F, user);}
2160:   else       {DMPlexComputeResidual_Internal(plex, cStart, cEnd, PETSC_MIN_REAL, X, NULL, F, user);}
2161:   DMDestroy(&plex);
2162:   return(0);
2163: }

2167: /*@
2168:   DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X

2170:   Input Parameters:
2171: + dm - The mesh
2172: - user - The user context

2174:   Output Parameter:
2175: . X  - Local solution

2177:   Level: developer

2179: .seealso: DMPlexComputeJacobianActionFEM()
2180: @*/
2181: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
2182: {
2183:   DM             plex;

2187:   DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2188:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
2189:   DMDestroy(&plex);
2190:   return(0);
2191: }

2195: PetscErrorCode DMPlexComputeJacobian_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user)
2196: {
2197:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
2198:   const char       *name  = "Jacobian";
2199:   DM                dmAux, plex;
2200:   DMLabel           depth;
2201:   Vec               A, cellgeom;
2202:   PetscDS           prob, probAux = NULL;
2203:   PetscQuadrature   quad;
2204:   PetscSection      section, globalSection, sectionAux;
2205:   PetscFECellGeom  *cgeom = NULL;
2206:   PetscScalar      *cgeomScal;
2207:   PetscScalar      *elemMat, *elemMatP, *elemMatD, *u, *u_t, *a = NULL;
2208:   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, c;
2209:   PetscInt          totDim, totDimBd, totDimAux, numBd, bd;
2210:   PetscBool         isShell, hasPrec, hasDyn;
2211:   PetscErrorCode    ierr;

2214:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
2215:   DMGetDimension(dm, &dim);
2216:   DMGetDefaultSection(dm, &section);
2217:   DMGetDefaultGlobalSection(dm, &globalSection);
2218:   DMGetDS(dm, &prob);
2219:   PetscDSGetTotalDimension(prob, &totDim);
2220:   PetscDSGetTotalBdDimension(prob, &totDimBd);
2221:   PetscDSHasJacobianPreconditioner(prob, &hasPrec);
2222:   PetscDSHasDynamicJacobian(prob, &hasDyn);
2223:   hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
2224:   PetscSectionGetNumFields(section, &Nf);
2225:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2226:   numCells = cEnd - cStart;
2227:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2228:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2229:   if (dmAux) {
2230:     DMConvert(dmAux, DMPLEX, &plex);
2231:     DMGetDefaultSection(plex, &sectionAux);
2232:     DMGetDS(dmAux, &probAux);
2233:     PetscDSGetTotalDimension(probAux, &totDimAux);
2234:   }
2235:   MatZeroEntries(JacP);
2236:   PetscMalloc5(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat,hasPrec ? numCells*totDim*totDim : 0, &elemMatP,hasDyn ? numCells*totDim*totDim : 0, &elemMatD);
2237:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2238:   DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2239:   VecGetArray(cellgeom, &cgeomScal);
2240:   if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2241:     DM dmCell;

2243:     VecGetDM(cellgeom,&dmCell);
2244:     PetscMalloc1(cEnd-cStart,&cgeom);
2245:     for (c = 0; c < cEnd - cStart; c++) {
2246:       PetscScalar *thisgeom;

2248:       DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2249:       cgeom[c] = *((PetscFECellGeom *) thisgeom);
2250:     }
2251:   }
2252:   else {
2253:     cgeom = (PetscFECellGeom *) cgeomScal;
2254:   }
2255:   for (c = cStart; c < cEnd; ++c) {
2256:     PetscScalar *x = NULL,  *x_t = NULL;
2257:     PetscInt     i;

2259:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2260:     for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2261:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2262:     if (X_t) {
2263:       DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2264:       for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2265:       DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2266:     }
2267:     if (dmAux) {
2268:       DMPlexVecGetClosure(plex, sectionAux, A, c, NULL, &x);
2269:       for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2270:       DMPlexVecRestoreClosure(plex, sectionAux, A, c, NULL, &x);
2271:     }
2272:   }
2273:   PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));
2274:   if (hasPrec) {PetscMemzero(elemMatP, numCells*totDim*totDim * sizeof(PetscScalar));}
2275:   if (hasDyn)  {PetscMemzero(elemMatD, numCells*totDim*totDim * sizeof(PetscScalar));}
2276:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
2277:     PetscFE  fe;
2278:     PetscInt numQuadPoints, Nb;
2279:     /* Conforming batches */
2280:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2281:     /* Remainder */
2282:     PetscInt Nr, offset;

2284:     PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2285:     PetscFEGetQuadrature(fe, &quad);
2286:     PetscFEGetDimension(fe, &Nb);
2287:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2288:     PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2289:     blockSize = Nb*numQuadPoints;
2290:     batchSize = numBlocks * blockSize;
2291:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2292:     numChunks = numCells / (numBatches*batchSize);
2293:     Ne        = numChunks*numBatches*batchSize;
2294:     Nr        = numCells % (numBatches*batchSize);
2295:     offset    = numCells - Nr;
2296:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2297:       PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMat);
2298:       PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMat[offset*totDim*totDim]);
2299:       if (hasPrec) {
2300:         PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMatP);
2301:         PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMatP[offset*totDim*totDim]);
2302:       }
2303:       if (hasDyn) {
2304:         PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMatD);
2305:         PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMatD[offset*totDim*totDim]);
2306:       }
2307:     }
2308:   }
2309:   if (hasDyn) {
2310:     for (c = 0; c < (cEnd - cStart)*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
2311:   }
2312:   for (c = cStart; c < cEnd; ++c) {
2313:     if (hasPrec) {
2314:       if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2315:       DMPlexMatSetClosure(dm, section, globalSection, Jac, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2316:       if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
2317:       DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMatP[(c-cStart)*totDim*totDim], ADD_VALUES);
2318:     } else {
2319:       if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2320:       DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2321:     }
2322:   }
2323:   if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2324:     PetscFree(cgeom);
2325:   }
2326:   else {
2327:     cgeom = NULL;
2328:   }
2329:   VecRestoreArray(cellgeom, &cgeomScal);
2330:   PetscFree5(u,u_t,elemMat,elemMatP,elemMatD);
2331:   if (dmAux) {
2332:     PetscFree(a);
2333:     DMDestroy(&plex);
2334:   }
2335:   DMPlexGetDepthLabel(dm, &depth);
2336:   DMGetNumBoundary(dm, &numBd);
2337:   DMPlexGetDepthLabel(dm, &depth);
2338:   DMGetNumBoundary(dm, &numBd);
2339:   for (bd = 0; bd < numBd; ++bd) {
2340:     const char     *bdLabel;
2341:     DMLabel         label;
2342:     IS              pointIS;
2343:     const PetscInt *points;
2344:     const PetscInt *values;
2345:     PetscInt        field, numValues, v, numPoints, p, dep, numFaces;
2346:     PetscBool       isEssential;
2347:     PetscObject     obj;
2348:     PetscClassId    id;

2350:     DMGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
2351:     DMGetField(dm, field, &obj);
2352:     PetscObjectGetClassId(obj, &id);
2353:     if ((id != PETSCFE_CLASSID) || isEssential) continue;
2354:     DMGetLabel(dm, bdLabel, &label);
2355:     for (v = 0; v < numValues; ++v) {
2356:       DMLabelGetStratumSize(label, values[v], &numPoints);
2357:       DMLabelGetStratumIS(label, values[v], &pointIS);
2358:       if (!pointIS) continue; /* No points with that id on this process */
2359:       ISGetIndices(pointIS, &points);
2360:       for (p = 0, numFaces = 0; p < numPoints; ++p) {
2361:         DMLabelGetValue(depth, points[p], &dep);
2362:         if (dep == dim-1) ++numFaces;
2363:       }
2364:       PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd*totDimBd,&elemMat);
2365:       if (X_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
2366:       for (p = 0, f = 0; p < numPoints; ++p) {
2367:         const PetscInt point = points[p];
2368:         PetscScalar   *x     = NULL;
2369:         PetscInt       i;

2371:         DMLabelGetValue(depth, points[p], &dep);
2372:         if (dep != dim-1) continue;
2373:         DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
2374:         DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
2375:         if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
2376:         DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
2377:         for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
2378:         DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
2379:         if (X_t) {
2380:           DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);
2381:           for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
2382:           DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);
2383:         }
2384:         ++f;
2385:       }
2386:       PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));
2387:       for (fieldI = 0; fieldI < Nf; ++fieldI) {
2388:         PetscFE  fe;
2389:         PetscInt numQuadPoints, Nb;
2390:         /* Conforming batches */
2391:         PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2392:         /* Remainder */
2393:         PetscInt Nr, offset;

2395:         PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
2396:         PetscFEGetQuadrature(fe, &quad);
2397:         PetscFEGetDimension(fe, &Nb);
2398:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2399:         PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2400:         blockSize = Nb*numQuadPoints;
2401:         batchSize = numBlocks * blockSize;
2402:          PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2403:         numChunks = numFaces / (numBatches*batchSize);
2404:         Ne        = numChunks*numBatches*batchSize;
2405:         Nr        = numFaces % (numBatches*batchSize);
2406:         offset    = numFaces - Nr;
2407:         for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2408:           PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, NULL, NULL, elemMat);
2409:           PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);
2410:         }
2411:       }
2412:       for (p = 0, f = 0; p < numPoints; ++p) {
2413:         const PetscInt point = points[p];

2415:         DMLabelGetValue(depth, point, &dep);
2416:         if (dep != dim-1) continue;
2417:         if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);}
2418:         DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);
2419:         ++f;
2420:       }
2421:       ISRestoreIndices(pointIS, &points);
2422:       ISDestroy(&pointIS);
2423:       PetscFree3(u,cgeom,elemMat);
2424:       if (X_t) {PetscFree(u_t);}
2425:     }
2426:   }
2427:   if (hasPrec) {
2428:     MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
2429:     MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
2430:   }
2431:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
2432:   MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
2433:   if (mesh->printFEM) {
2434:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
2435:     MatChop(JacP, 1.0e-10);
2436:     MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
2437:   }
2438:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2439:   PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
2440:   if (isShell) {
2441:     JacActionCtx *jctx;

2443:     MatShellGetContext(Jac, &jctx);
2444:     VecCopy(X, jctx->u);
2445:   }
2446:   return(0);
2447: }

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

2454:   Input Parameters:
2455: + dm - The mesh
2456: . X  - Local input vector
2457: - user - The user context

2459:   Output Parameter:
2460: . Jac  - Jacobian matrix

2462:   Note:
2463:   The first member of the user context must be an FEMContext.

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

2468:   Level: developer

2470: .seealso: FormFunctionLocal()
2471: @*/
2472: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
2473: {
2474:   PetscInt       cStart, cEnd, cEndInterior;
2475:   DM             plex;

2479:   DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2480:   DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2481:   DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2482:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2483:   DMPlexComputeJacobian_Internal(plex, cStart, cEnd, 0.0, 0.0, X, NULL, Jac, JacP, user);
2484:   DMDestroy(&plex);
2485:   return(0);
2486: }

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

2493:   Input Parameters:
2494: + dm - The DM object
2495: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see DMAddBoundary())
2496: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
2497: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())

2499:   Level: developer
2500: @*/
2501: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
2502: {

2506:   DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);
2507:   DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);
2508:   DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);
2509:   return(0);
2510: }

2514: PetscErrorCode DMSNESCheckFromOptions_Internal(SNES snes, DM dm, Vec u, Vec sol, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx), void **ctxs)
2515: {
2516:   Mat            J, M;
2517:   Vec            r, b;
2518:   MatNullSpace   nullSpace;
2519:   PetscReal     *error, res = 0.0;
2520:   PetscInt       numFields;

2524:   VecDuplicate(u, &r);
2525:   DMCreateMatrix(dm, &J);
2526:   M    = J;
2527:   /* TODO Null space for J */
2528:   /* Check discretization error */
2529:   DMGetNumFields(dm, &numFields);
2530:   PetscMalloc1(PetscMax(1, numFields), &error);
2531:   if (numFields > 1) {
2532:     PetscInt f;

2534:     DMComputeL2FieldDiff(dm, 0.0, exactFuncs, ctxs, u, error);
2535:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: [");
2536:     for (f = 0; f < numFields; ++f) {
2537:       if (f) {PetscPrintf(PETSC_COMM_WORLD, ", ");}
2538:       if (error[f] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "%g", error[f]);}
2539:       else                     {PetscPrintf(PETSC_COMM_WORLD, "< 1.0e-11");}
2540:     }
2541:     PetscPrintf(PETSC_COMM_WORLD, "]\n");
2542:   } else {
2543:     DMComputeL2Diff(dm, 0.0, exactFuncs, ctxs, u, &error[0]);
2544:     if (error[0] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error[0]);}
2545:     else                     {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
2546:   }
2547:   PetscFree(error);
2548:   /* Check residual */
2549:   SNESComputeFunction(snes, u, r);
2550:   VecNorm(r, NORM_2, &res);
2551:   PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
2552:   VecChop(r, 1.0e-10);
2553:   PetscObjectSetName((PetscObject) r, "Initial Residual");
2554:   PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
2555:   VecViewFromOptions(r, NULL, "-vec_view");
2556:   /* Check Jacobian */
2557:   SNESComputeJacobian(snes, u, M, M);
2558:   MatGetNullSpace(J, &nullSpace);
2559:   if (nullSpace) {
2560:     PetscBool isNull;
2561:     MatNullSpaceTest(nullSpace, J, &isNull);
2562:     if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
2563:   }
2564:   VecDuplicate(u, &b);
2565:   VecSet(r, 0.0);
2566:   SNESComputeFunction(snes, r, b);
2567:   MatMult(M, u, r);
2568:   VecAXPY(r, 1.0, b);
2569:   VecDestroy(&b);
2570:   VecNorm(r, NORM_2, &res);
2571:   PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
2572:   VecChop(r, 1.0e-10);
2573:   PetscObjectSetName((PetscObject) r, "Au - b = Au + F(0)");
2574:   PetscObjectSetOptionsPrefix((PetscObject)r,"linear_res_");
2575:   VecViewFromOptions(r, NULL, "-vec_view");
2576:   VecDestroy(&r);
2577:   MatNullSpaceDestroy(&nullSpace);
2578:   MatDestroy(&J);
2579:   return(0);
2580: }

2584: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
2585: {
2586:   DM             dm;
2587:   Vec            sol;
2588:   PetscBool      check;

2592:   PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);
2593:   if (!check) return(0);
2594:   SNESGetDM(snes, &dm);
2595:   VecDuplicate(u, &sol);
2596:   SNESSetSolution(snes, sol);
2597:   DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
2598:   VecDestroy(&sol);
2599:   return(0);
2600: }