Actual source code: dmplexsnes.c

petsc-3.6.4 2016-04-12
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: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
 12: {

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

 19:   (*ctx)->comm   = comm;
 20:   (*ctx)->dim    = -1;
 21:   (*ctx)->nInput = 0;
 22:   (*ctx)->points = NULL;
 23:   (*ctx)->cells  = NULL;
 24:   (*ctx)->n      = -1;
 25:   (*ctx)->coords = NULL;
 26:   return(0);
 27: }

 31: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
 32: {
 34:   if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);
 35:   ctx->dim = dim;
 36:   return(0);
 37: }

 41: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
 42: {
 45:   *dim = ctx->dim;
 46:   return(0);
 47: }

 51: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
 52: {
 54:   if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);
 55:   ctx->dof = dof;
 56:   return(0);
 57: }

 61: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
 62: {
 65:   *dof = ctx->dof;
 66:   return(0);
 67: }

 71: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
 72: {

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

 80:   PetscMalloc1(n*ctx->dim, &ctx->points);
 81:   PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));
 82:   return(0);
 83: }

 87: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
 88: {
 89:   MPI_Comm       comm = ctx->comm;
 90:   PetscScalar   *a;
 91:   PetscInt       p, q, i;
 92:   PetscMPIInt    rank, size;
 94:   Vec            pointVec;
 95:   IS             cellIS;
 96:   PetscLayout    layout;
 97:   PetscReal      *globalPoints;
 98:   PetscScalar    *globalPointsScalar;
 99:   const PetscInt *ranges;
100:   PetscMPIInt    *counts, *displs;
101:   const PetscInt *foundCells;
102:   PetscMPIInt    *foundProcs, *globalProcs;
103:   PetscInt       n, N;

107:   MPI_Comm_size(comm, &size);
108:   MPI_Comm_rank(comm, &rank);
109:   if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
110:   /* Locate points */
111:   n = ctx->nInput;
112:   if (!redundantPoints) {
113:     PetscLayoutCreate(comm, &layout);
114:     PetscLayoutSetBlockSize(layout, 1);
115:     PetscLayoutSetLocalSize(layout, n);
116:     PetscLayoutSetUp(layout);
117:     PetscLayoutGetSize(layout, &N);
118:     /* Communicate all points to all processes */
119:     PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
120:     PetscLayoutGetRanges(layout, &ranges);
121:     for (p = 0; p < size; ++p) {
122:       counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
123:       displs[p] = ranges[p]*ctx->dim;
124:     }
125:     MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
126:   } else {
127:     N = n;
128:     globalPoints = ctx->points;
129:     counts = displs = NULL;
130:     layout = NULL;
131:   }
132: #if 0
133:   PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
134:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
135: #else
136: #if defined(PETSC_USE_COMPLEX)
137:   PetscMalloc1(N,&globalPointsScalar);
138:   for (i=0; i<N; i++) globalPointsScalar[i] = globalPoints[i];
139: #else
140:   globalPointsScalar = globalPoints;
141: #endif
142:   VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
143:   PetscMalloc2(N,&foundProcs,N,&globalProcs);
144:   DMLocatePoints(dm, pointVec, &cellIS);
145:   ISGetIndices(cellIS, &foundCells);
146: #endif
147:   for (p = 0; p < N; ++p) {
148:     if (foundCells[p] >= 0) foundProcs[p] = rank;
149:     else foundProcs[p] = size;
150:   }
151:   /* Let the lowest rank process own each point */
152:   MPI_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
153:   ctx->n = 0;
154:   for (p = 0; p < N; ++p) {
155:     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);
156:     else if (globalProcs[p] == rank) ctx->n++;
157:   }
158:   /* Create coordinates vector and array of owned cells */
159:   PetscMalloc1(ctx->n, &ctx->cells);
160:   VecCreate(comm, &ctx->coords);
161:   VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
162:   VecSetBlockSize(ctx->coords, ctx->dim);
163:   VecSetType(ctx->coords,VECSTANDARD);
164:   VecGetArray(ctx->coords, &a);
165:   for (p = 0, q = 0, i = 0; p < N; ++p) {
166:     if (globalProcs[p] == rank) {
167:       PetscInt d;

169:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
170:       ctx->cells[q++] = foundCells[p];
171:     }
172:   }
173:   VecRestoreArray(ctx->coords, &a);
174: #if 0
175:   PetscFree3(foundCells,foundProcs,globalProcs);
176: #else
177:   PetscFree2(foundProcs,globalProcs);
178:   ISRestoreIndices(cellIS, &foundCells);
179:   ISDestroy(&cellIS);
180:   VecDestroy(&pointVec);
181: #endif
182:   if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
183:   if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
184:   PetscLayoutDestroy(&layout);
185:   return(0);
186: }

190: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
191: {
194:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
195:   *coordinates = ctx->coords;
196:   return(0);
197: }

201: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
202: {

207:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
208:   VecCreate(ctx->comm, v);
209:   VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
210:   VecSetBlockSize(*v, ctx->dof);
211:   VecSetType(*v,VECSTANDARD);
212:   return(0);
213: }

217: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
218: {

223:   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
224:   VecDestroy(v);
225:   return(0);
226: }

230: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
231: {
232:   PetscReal      *v0, *J, *invJ, detJ;
233:   const PetscScalar *coords;
234:   PetscScalar    *a;
235:   PetscInt       p;

239:   PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
240:   VecGetArrayRead(ctx->coords, &coords);
241:   VecGetArray(v, &a);
242:   for (p = 0; p < ctx->n; ++p) {
243:     PetscInt     c = ctx->cells[p];
244:     PetscScalar *x = NULL;
245:     PetscReal    xi[4];
246:     PetscInt     d, f, comp;

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

253:     for (d = 0; d < ctx->dim; ++d) {
254:       xi[d] = 0.0;
255:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
256:       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];
257:     }
258:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
259:   }
260:   VecRestoreArray(v, &a);
261:   VecRestoreArrayRead(ctx->coords, &coords);
262:   PetscFree3(v0, J, invJ);
263:   return(0);
264: }

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

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

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

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

307: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
308: {
309:   const PetscScalar *vertices = (const PetscScalar*) ctx;
310:   const PetscScalar x0        = vertices[0];
311:   const PetscScalar y0        = vertices[1];
312:   const PetscScalar x1        = vertices[2];
313:   const PetscScalar y1        = vertices[3];
314:   const PetscScalar x2        = vertices[4];
315:   const PetscScalar y2        = vertices[5];
316:   const PetscScalar x3        = vertices[6];
317:   const PetscScalar y3        = vertices[7];
318:   const PetscScalar f_1       = x1 - x0;
319:   const PetscScalar g_1       = y1 - y0;
320:   const PetscScalar f_3       = x3 - x0;
321:   const PetscScalar g_3       = y3 - y0;
322:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
323:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
324:   const PetscScalar *ref;
325:   PetscScalar       *real;
326:   PetscErrorCode    ierr;

329:   VecGetArrayRead(Xref,  &ref);
330:   VecGetArray(Xreal, &real);
331:   {
332:     const PetscScalar p0 = ref[0];
333:     const PetscScalar p1 = ref[1];

335:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
336:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
337:   }
338:   PetscLogFlops(28);
339:   VecRestoreArrayRead(Xref,  &ref);
340:   VecRestoreArray(Xreal, &real);
341:   return(0);
342: }

344: #include <petsc/private/dmimpl.h>
347: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
348: {
349:   const PetscScalar *vertices = (const PetscScalar*) ctx;
350:   const PetscScalar x0        = vertices[0];
351:   const PetscScalar y0        = vertices[1];
352:   const PetscScalar x1        = vertices[2];
353:   const PetscScalar y1        = vertices[3];
354:   const PetscScalar x2        = vertices[4];
355:   const PetscScalar y2        = vertices[5];
356:   const PetscScalar x3        = vertices[6];
357:   const PetscScalar y3        = vertices[7];
358:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
359:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
360:   const PetscScalar *ref;
361:   PetscErrorCode    ierr;

364:   VecGetArrayRead(Xref,  &ref);
365:   {
366:     const PetscScalar x       = ref[0];
367:     const PetscScalar y       = ref[1];
368:     const PetscInt    rows[2] = {0, 1};
369:     PetscScalar       values[4];

371:     values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
372:     values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
373:     MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
374:   }
375:   PetscLogFlops(30);
376:   VecRestoreArrayRead(Xref,  &ref);
377:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
378:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
379:   return(0);
380: }

384: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
385: {
386:   DM             dmCoord;
387:   SNES           snes;
388:   KSP            ksp;
389:   PC             pc;
390:   Vec            coordsLocal, r, ref, real;
391:   Mat            J;
392:   const PetscScalar *coords;
393:   PetscScalar    *a;
394:   PetscInt       p;

398:   DMGetCoordinatesLocal(dm, &coordsLocal);
399:   DMGetCoordinateDM(dm, &dmCoord);
400:   SNESCreate(PETSC_COMM_SELF, &snes);
401:   SNESSetOptionsPrefix(snes, "quad_interp_");
402:   VecCreate(PETSC_COMM_SELF, &r);
403:   VecSetSizes(r, 2, 2);
404:   VecSetType(r,dm->vectype);
405:   VecDuplicate(r, &ref);
406:   VecDuplicate(r, &real);
407:   MatCreate(PETSC_COMM_SELF, &J);
408:   MatSetSizes(J, 2, 2, 2, 2);
409:   MatSetType(J, MATSEQDENSE);
410:   MatSetUp(J);
411:   SNESSetFunction(snes, r, QuadMap_Private, NULL);
412:   SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
413:   SNESGetKSP(snes, &ksp);
414:   KSPGetPC(ksp, &pc);
415:   PCSetType(pc, PCLU);
416:   SNESSetFromOptions(snes);

418:   VecGetArrayRead(ctx->coords, &coords);
419:   VecGetArray(v, &a);
420:   for (p = 0; p < ctx->n; ++p) {
421:     PetscScalar *x = NULL, *vertices = NULL;
422:     PetscScalar *xi;
423:     PetscReal    xir[2];
424:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

426:     /* Can make this do all points at once */
427:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
428:     if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);
429:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
430:     if (4*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 4*ctx->dof);
431:     SNESSetFunction(snes, NULL, NULL, (void*) vertices);
432:     SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
433:     VecGetArray(real, &xi);
434:     xi[0]  = coords[p*ctx->dim+0];
435:     xi[1]  = coords[p*ctx->dim+1];
436:     VecRestoreArray(real, &xi);
437:     SNESSolve(snes, real, ref);
438:     VecGetArray(ref, &xi);
439:     xir[0] = PetscRealPart(xi[0]);
440:     xir[1] = PetscRealPart(xi[1]);
441:     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];

443:     VecRestoreArray(ref, &xi);
444:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
445:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
446:   }
447:   VecRestoreArray(v, &a);
448:   VecRestoreArrayRead(ctx->coords, &coords);

450:   SNESDestroy(&snes);
451:   VecDestroy(&r);
452:   VecDestroy(&ref);
453:   VecDestroy(&real);
454:   MatDestroy(&J);
455:   return(0);
456: }

460: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
461: {
462:   const PetscScalar *vertices = (const PetscScalar*) ctx;
463:   const PetscScalar x0        = vertices[0];
464:   const PetscScalar y0        = vertices[1];
465:   const PetscScalar z0        = vertices[2];
466:   const PetscScalar x1        = vertices[9];
467:   const PetscScalar y1        = vertices[10];
468:   const PetscScalar z1        = vertices[11];
469:   const PetscScalar x2        = vertices[6];
470:   const PetscScalar y2        = vertices[7];
471:   const PetscScalar z2        = vertices[8];
472:   const PetscScalar x3        = vertices[3];
473:   const PetscScalar y3        = vertices[4];
474:   const PetscScalar z3        = vertices[5];
475:   const PetscScalar x4        = vertices[12];
476:   const PetscScalar y4        = vertices[13];
477:   const PetscScalar z4        = vertices[14];
478:   const PetscScalar x5        = vertices[15];
479:   const PetscScalar y5        = vertices[16];
480:   const PetscScalar z5        = vertices[17];
481:   const PetscScalar x6        = vertices[18];
482:   const PetscScalar y6        = vertices[19];
483:   const PetscScalar z6        = vertices[20];
484:   const PetscScalar x7        = vertices[21];
485:   const PetscScalar y7        = vertices[22];
486:   const PetscScalar z7        = vertices[23];
487:   const PetscScalar f_1       = x1 - x0;
488:   const PetscScalar g_1       = y1 - y0;
489:   const PetscScalar h_1       = z1 - z0;
490:   const PetscScalar f_3       = x3 - x0;
491:   const PetscScalar g_3       = y3 - y0;
492:   const PetscScalar h_3       = z3 - z0;
493:   const PetscScalar f_4       = x4 - x0;
494:   const PetscScalar g_4       = y4 - y0;
495:   const PetscScalar h_4       = z4 - z0;
496:   const PetscScalar f_01      = x2 - x1 - x3 + x0;
497:   const PetscScalar g_01      = y2 - y1 - y3 + y0;
498:   const PetscScalar h_01      = z2 - z1 - z3 + z0;
499:   const PetscScalar f_12      = x7 - x3 - x4 + x0;
500:   const PetscScalar g_12      = y7 - y3 - y4 + y0;
501:   const PetscScalar h_12      = z7 - z3 - z4 + z0;
502:   const PetscScalar f_02      = x5 - x1 - x4 + x0;
503:   const PetscScalar g_02      = y5 - y1 - y4 + y0;
504:   const PetscScalar h_02      = z5 - z1 - z4 + z0;
505:   const PetscScalar f_012     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
506:   const PetscScalar g_012     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
507:   const PetscScalar h_012     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
508:   const PetscScalar *ref;
509:   PetscScalar       *real;
510:   PetscErrorCode    ierr;

513:   VecGetArrayRead(Xref,  &ref);
514:   VecGetArray(Xreal, &real);
515:   {
516:     const PetscScalar p0 = ref[0];
517:     const PetscScalar p1 = ref[1];
518:     const PetscScalar p2 = ref[2];

520:     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;
521:     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;
522:     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;
523:   }
524:   PetscLogFlops(114);
525:   VecRestoreArrayRead(Xref,  &ref);
526:   VecRestoreArray(Xreal, &real);
527:   return(0);
528: }

532: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
533: {
534:   const PetscScalar *vertices = (const PetscScalar*) ctx;
535:   const PetscScalar x0        = vertices[0];
536:   const PetscScalar y0        = vertices[1];
537:   const PetscScalar z0        = vertices[2];
538:   const PetscScalar x1        = vertices[9];
539:   const PetscScalar y1        = vertices[10];
540:   const PetscScalar z1        = vertices[11];
541:   const PetscScalar x2        = vertices[6];
542:   const PetscScalar y2        = vertices[7];
543:   const PetscScalar z2        = vertices[8];
544:   const PetscScalar x3        = vertices[3];
545:   const PetscScalar y3        = vertices[4];
546:   const PetscScalar z3        = vertices[5];
547:   const PetscScalar x4        = vertices[12];
548:   const PetscScalar y4        = vertices[13];
549:   const PetscScalar z4        = vertices[14];
550:   const PetscScalar x5        = vertices[15];
551:   const PetscScalar y5        = vertices[16];
552:   const PetscScalar z5        = vertices[17];
553:   const PetscScalar x6        = vertices[18];
554:   const PetscScalar y6        = vertices[19];
555:   const PetscScalar z6        = vertices[20];
556:   const PetscScalar x7        = vertices[21];
557:   const PetscScalar y7        = vertices[22];
558:   const PetscScalar z7        = vertices[23];
559:   const PetscScalar f_xy      = x2 - x1 - x3 + x0;
560:   const PetscScalar g_xy      = y2 - y1 - y3 + y0;
561:   const PetscScalar h_xy      = z2 - z1 - z3 + z0;
562:   const PetscScalar f_yz      = x7 - x3 - x4 + x0;
563:   const PetscScalar g_yz      = y7 - y3 - y4 + y0;
564:   const PetscScalar h_yz      = z7 - z3 - z4 + z0;
565:   const PetscScalar f_xz      = x5 - x1 - x4 + x0;
566:   const PetscScalar g_xz      = y5 - y1 - y4 + y0;
567:   const PetscScalar h_xz      = z5 - z1 - z4 + z0;
568:   const PetscScalar f_xyz     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
569:   const PetscScalar g_xyz     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
570:   const PetscScalar h_xyz     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
571:   const PetscScalar *ref;
572:   PetscErrorCode    ierr;

575:   VecGetArrayRead(Xref,  &ref);
576:   {
577:     const PetscScalar x       = ref[0];
578:     const PetscScalar y       = ref[1];
579:     const PetscScalar z       = ref[2];
580:     const PetscInt    rows[3] = {0, 1, 2};
581:     PetscScalar       values[9];

583:     values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
584:     values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
585:     values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
586:     values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
587:     values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
588:     values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
589:     values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
590:     values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
591:     values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;

593:     MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
594:   }
595:   PetscLogFlops(152);
596:   VecRestoreArrayRead(Xref,  &ref);
597:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
598:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
599:   return(0);
600: }

604: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
605: {
606:   DM             dmCoord;
607:   SNES           snes;
608:   KSP            ksp;
609:   PC             pc;
610:   Vec            coordsLocal, r, ref, real;
611:   Mat            J;
612:   const PetscScalar *coords;
613:   PetscScalar    *a;
614:   PetscInt       p;

618:   DMGetCoordinatesLocal(dm, &coordsLocal);
619:   DMGetCoordinateDM(dm, &dmCoord);
620:   SNESCreate(PETSC_COMM_SELF, &snes);
621:   SNESSetOptionsPrefix(snes, "hex_interp_");
622:   VecCreate(PETSC_COMM_SELF, &r);
623:   VecSetSizes(r, 3, 3);
624:   VecSetType(r,dm->vectype);
625:   VecDuplicate(r, &ref);
626:   VecDuplicate(r, &real);
627:   MatCreate(PETSC_COMM_SELF, &J);
628:   MatSetSizes(J, 3, 3, 3, 3);
629:   MatSetType(J, MATSEQDENSE);
630:   MatSetUp(J);
631:   SNESSetFunction(snes, r, HexMap_Private, NULL);
632:   SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
633:   SNESGetKSP(snes, &ksp);
634:   KSPGetPC(ksp, &pc);
635:   PCSetType(pc, PCLU);
636:   SNESSetFromOptions(snes);

638:   VecGetArrayRead(ctx->coords, &coords);
639:   VecGetArray(v, &a);
640:   for (p = 0; p < ctx->n; ++p) {
641:     PetscScalar *x = NULL, *vertices = NULL;
642:     PetscScalar *xi;
643:     PetscReal    xir[3];
644:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

646:     /* Can make this do all points at once */
647:     DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
648:     if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);
649:     DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
650:     if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);
651:     SNESSetFunction(snes, NULL, NULL, (void*) vertices);
652:     SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
653:     VecGetArray(real, &xi);
654:     xi[0]  = coords[p*ctx->dim+0];
655:     xi[1]  = coords[p*ctx->dim+1];
656:     xi[2]  = coords[p*ctx->dim+2];
657:     VecRestoreArray(real, &xi);
658:     SNESSolve(snes, real, ref);
659:     VecGetArray(ref, &xi);
660:     xir[0] = PetscRealPart(xi[0]);
661:     xir[1] = PetscRealPart(xi[1]);
662:     xir[2] = PetscRealPart(xi[2]);
663:     for (comp = 0; comp < ctx->dof; ++comp) {
664:       a[p*ctx->dof+comp] =
665:         x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
666:         x[3*ctx->dof+comp]*    xir[0]*(1-xir[1])*(1-xir[2]) +
667:         x[2*ctx->dof+comp]*    xir[0]*    xir[1]*(1-xir[2]) +
668:         x[1*ctx->dof+comp]*(1-xir[0])*    xir[1]*(1-xir[2]) +
669:         x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*   xir[2] +
670:         x[5*ctx->dof+comp]*    xir[0]*(1-xir[1])*   xir[2] +
671:         x[6*ctx->dof+comp]*    xir[0]*    xir[1]*   xir[2] +
672:         x[7*ctx->dof+comp]*(1-xir[0])*    xir[1]*   xir[2];
673:     }
674:     VecRestoreArray(ref, &xi);
675:     DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
676:     DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
677:   }
678:   VecRestoreArray(v, &a);
679:   VecRestoreArrayRead(ctx->coords, &coords);

681:   SNESDestroy(&snes);
682:   VecDestroy(&r);
683:   VecDestroy(&ref);
684:   VecDestroy(&real);
685:   MatDestroy(&J);
686:   return(0);
687: }

691: /*
692:   Input Parameters:
693: + ctx - The DMInterpolationInfo context
694: . dm  - The DM
695: - x   - The local vector containing the field to be interpolated

697:   Output Parameters:
698: . v   - The vector containing the interpolated values
699: */
700: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
701: {
702:   PetscInt       dim, coneSize, n;

709:   VecGetLocalSize(v, &n);
710:   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);
711:   if (n) {
712:     DMGetDimension(dm, &dim);
713:     DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);
714:     if (dim == 2) {
715:       if (coneSize == 3) {
716:         DMInterpolate_Triangle_Private(ctx, dm, x, v);
717:       } else if (coneSize == 4) {
718:         DMInterpolate_Quad_Private(ctx, dm, x, v);
719:       } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
720:     } else if (dim == 3) {
721:       if (coneSize == 4) {
722:         DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
723:       } else {
724:         DMInterpolate_Hex_Private(ctx, dm, x, v);
725:       }
726:     } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
727:   }
728:   return(0);
729: }

733: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
734: {

739:   VecDestroy(&(*ctx)->coords);
740:   PetscFree((*ctx)->points);
741:   PetscFree((*ctx)->cells);
742:   PetscFree(*ctx);
743:   *ctx = NULL;
744:   return(0);
745: }

749: /*@C
750:   SNESMonitorFields - Monitors the residual for each field separately

752:   Collective on SNES

754:   Input Parameters:
755: + snes   - the SNES context
756: . its    - iteration number
757: . fgnorm - 2-norm of residual
758: - dummy  - unused context

760:   Notes:
761:   This routine prints the residual norm at each iteration.

763:   Level: intermediate

765: .keywords: SNES, nonlinear, default, monitor, norm
766: .seealso: SNESMonitorSet(), SNESMonitorDefault()
767: @*/
768: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
769: {
770:   PetscViewer        viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) snes));
771:   Vec                res;
772:   DM                 dm;
773:   PetscSection       s;
774:   const PetscScalar *r;
775:   PetscReal         *lnorms, *norms;
776:   PetscInt           numFields, f, pStart, pEnd, p;
777:   PetscErrorCode     ierr;

780:   SNESGetFunction(snes, &res, 0, 0);
781:   SNESGetDM(snes, &dm);
782:   DMGetDefaultSection(dm, &s);
783:   PetscSectionGetNumFields(s, &numFields);
784:   PetscSectionGetChart(s, &pStart, &pEnd);
785:   PetscCalloc2(numFields, &lnorms, numFields, &norms);
786:   VecGetArrayRead(res, &r);
787:   for (p = pStart; p < pEnd; ++p) {
788:     for (f = 0; f < numFields; ++f) {
789:       PetscInt fdof, foff, d;

791:       PetscSectionGetFieldDof(s, p, f, &fdof);
792:       PetscSectionGetFieldOffset(s, p, f, &foff);
793:       for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
794:     }
795:   }
796:   VecRestoreArrayRead(res, &r);
797:   MPI_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
798:   PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
799:   PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
800:   for (f = 0; f < numFields; ++f) {
801:     if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
802:     PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
803:   }
804:   PetscViewerASCIIPrintf(viewer, "]\n");
805:   PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
806:   PetscFree2(lnorms, norms);
807:   return(0);
808: }

810: /********************* Residual Computation **************************/

814: /*@
815:   DMPlexSNESGetGeometryFEM - Return precomputed geometric data

817:   Input Parameter:
818: . dm - The DM

820:   Output Parameters:
821: . cellgeom - The values precomputed from cell geometry

823:   Level: developer

825: .seealso: DMPlexSNESSetFunctionLocal()
826: @*/
827: PetscErrorCode DMPlexSNESGetGeometryFEM(DM dm, Vec *cellgeom)
828: {
829:   DMSNES         dmsnes;
830:   PetscObject    obj;

835:   DMGetDMSNES(dm, &dmsnes);
836:   PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", &obj);
837:   if (!obj) {
838:     Vec cellgeom;

840:     DMPlexComputeGeometryFEM(dm, &cellgeom);
841:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", (PetscObject) cellgeom);
842:     VecDestroy(&cellgeom);
843:   }
845:   return(0);
846: }

850: /*@
851:   DMPlexSNESGetGeometryFVM - Return precomputed geometric data

853:   Input Parameter:
854: . dm - The DM

856:   Output Parameters:
857: + facegeom - The values precomputed from face geometry
858: . cellgeom - The values precomputed from cell geometry
859: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell

861:   Level: developer

863: .seealso: DMPlexTSSetRHSFunctionLocal()
864: @*/
865: PetscErrorCode DMPlexSNESGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
866: {
867:   DMSNES         dmsnes;
868:   PetscObject    obj;

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

878:     DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);
879:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_facegeom_fvm", (PetscObject) facegeom);
880:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fvm", (PetscObject) cellgeom);
881:     VecDestroy(&facegeom);
882:     VecDestroy(&cellgeom);
883:   }
886:   if (minRadius) {DMPlexGetMinRadius(dm, minRadius);}
887:   return(0);
888: }

892: /*@
893:   DMPlexSNESGetGradientDM - Return gradient data layout

895:   Input Parameters:
896: + dm - The DM
897: - fv - The PetscFV

899:   Output Parameter:
900: . dmGrad - The layout for gradient values

902:   Level: developer

904: .seealso: DMPlexSNESGetGeometryFVM()
905: @*/
906: PetscErrorCode DMPlexSNESGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
907: {
908:   DMSNES         dmsnes;
909:   PetscObject    obj;
910:   PetscBool      computeGradients;

917:   PetscFVGetComputeGradients(fv, &computeGradients);
918:   if (!computeGradients) {*dmGrad = NULL; return(0);}
919:   DMGetDMSNES(dm, &dmsnes);
920:   PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", &obj);
921:   if (!obj) {
922:     DM  dmGrad;
923:     Vec faceGeometry, cellGeometry;

925:     DMPlexSNESGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);
926:     DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);
927:     PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject) dmGrad);
928:     DMDestroy(&dmGrad);
929:   }
930:   PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject *) dmGrad);
931:   return(0);
932: }

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

939:   Input Parameters:
940: + dm     - The DM
941: . cStart - The first cell to include
942: . cEnd   - The first cell to exclude
943: . locX   - A local vector with the solution fields
944: . locX_t - A local vector with solution field time derivatives, or NULL
945: - locA   - A local vector with auxiliary fields, or NULL

947:   Output Parameters:
948: + u   - The field coefficients
949: . u_t - The fields derivative coefficients
950: - a   - The auxiliary field coefficients

952:   Level: developer

954: .seealso: DMPlexGetFaceFields()
955: @*/
956: PetscErrorCode DMPlexGetCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
957: {
958:   DM             dmAux;
959:   PetscSection   section, sectionAux;
960:   PetscDS        prob;
961:   PetscInt       numCells = cEnd - cStart, totDim, totDimAux, c;

972:   DMGetDefaultSection(dm, &section);
973:   DMGetDS(dm, &prob);
974:   PetscDSGetTotalDimension(prob, &totDim);
975:   if (locA) {
976:     PetscDS probAux;

978:     VecGetDM(locA, &dmAux);
979:     DMGetDefaultSection(dmAux, &sectionAux);
980:     DMGetDS(dmAux, &probAux);
981:     PetscDSGetTotalDimension(probAux, &totDimAux);
982:   }
983:   DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u);
984:   if (locX_t) {DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u_t);} else {*u_t = NULL;}
985:   if (locA)   {DMGetWorkArray(dm, numCells*totDimAux, PETSC_SCALAR, a);} else {*a = NULL;}
986:   for (c = cStart; c < cEnd; ++c) {
987:     PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
988:     PetscInt     i;

990:     DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
991:     for (i = 0; i < totDim; ++i) ul[(c-cStart)*totDim+i] = x[i];
992:     DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
993:     if (locX_t) {
994:       DMPlexVecGetClosure(dm, section, locX_t, c, NULL, &x_t);
995:       for (i = 0; i < totDim; ++i) ul_t[(c-cStart)*totDim+i] = x_t[i];
996:       DMPlexVecRestoreClosure(dm, section, locX_t, c, NULL, &x_t);
997:     }
998:     if (locA) {
999:       DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1000:       for (i = 0; i < totDimAux; ++i) al[(c-cStart)*totDimAux+i] = x[i];
1001:       DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1002:     }
1003:   }
1004:   return(0);
1005: }

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

1012:   Input Parameters:
1013: + dm     - The DM
1014: . cStart - The first cell to include
1015: . cEnd   - The first cell to exclude
1016: . locX   - A local vector with the solution fields
1017: . locX_t - A local vector with solution field time derivatives, or NULL
1018: - locA   - A local vector with auxiliary fields, or NULL

1020:   Output Parameters:
1021: + u   - The field coefficients
1022: . u_t - The fields derivative coefficients
1023: - a   - The auxiliary field coefficients

1025:   Level: developer

1027: .seealso: DMPlexGetFaceFields()
1028: @*/
1029: PetscErrorCode DMPlexRestoreCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
1030: {

1034:   DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u);
1035:   if (*u_t) {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u_t);}
1036:   if (*a)   {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, a);}
1037:   return(0);
1038: }

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

1045:   Input Parameters:
1046: + dm     - The DM
1047: . fStart - The first face to include
1048: . fEnd   - The first face to exclude
1049: . locX   - A local vector with the solution fields
1050: . locX_t - A local vector with solution field time derivatives, or NULL
1051: . faceGeometry - A local vector with face geometry
1052: . cellGeometry - A local vector with cell geometry
1053: - locaGrad - A local vector with field gradients, or NULL

1055:   Output Parameters:
1056: + uL - The field values at the left side of the face
1057: - uR - The field values at the right side of the face

1059:   Level: developer

1061: .seealso: DMPlexGetCellFields()
1062: @*/
1063: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1064: {
1065:   DM                 dmFace, dmCell, dmGrad = NULL;
1066:   PetscSection       section;
1067:   PetscDS            prob;
1068:   DMLabel            ghostLabel;
1069:   const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
1070:   PetscBool         *isFE;
1071:   PetscInt           dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
1072:   PetscErrorCode     ierr;

1083:   DMGetDimension(dm, &dim);
1084:   DMGetDS(dm, &prob);
1085:   DMGetDefaultSection(dm, &section);
1086:   PetscDSGetNumFields(prob, &Nf);
1087:   PetscDSGetTotalComponents(prob, &Nc);
1088:   PetscMalloc1(Nf, &isFE);
1089:   for (f = 0; f < Nf; ++f) {
1090:     PetscObject  obj;
1091:     PetscClassId id;

1093:     DMGetField(dm, f, &obj);
1094:     PetscObjectGetClassId(obj, &id);
1095:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
1096:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
1097:     else                            {isFE[f] = PETSC_FALSE;}
1098:   }
1099:   DMPlexGetLabel(dm, "ghost", &ghostLabel);
1100:   VecGetArrayRead(locX, &x);
1101:   VecGetDM(faceGeometry, &dmFace);
1102:   VecGetArrayRead(faceGeometry, &facegeom);
1103:   VecGetDM(cellGeometry, &dmCell);
1104:   VecGetArrayRead(cellGeometry, &cellgeom);
1105:   if (locGrad) {
1106:     VecGetDM(locGrad, &dmGrad);
1107:     VecGetArrayRead(locGrad, &lgrad);
1108:   }
1109:   DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uL);
1110:   DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uR);
1111:   /* Right now just eat the extra work for FE (could make a cell loop) */
1112:   for (face = fStart, iface = 0; face < fEnd; ++face) {
1113:     const PetscInt        *cells;
1114:     const PetscFVFaceGeom *fg;
1115:     const PetscFVCellGeom *cgL, *cgR;
1116:     const PetscScalar     *xL, *xR, *gL, *gR;
1117:     PetscScalar           *uLl = *uL, *uRl = *uR;
1118:     PetscInt               ghost;

1120:     DMLabelGetValue(ghostLabel, face, &ghost);
1121:     if (ghost >= 0) continue;
1122:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1123:     DMPlexGetSupport(dm, face, &cells);
1124:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1125:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1126:     for (f = 0; f < Nf; ++f) {
1127:       PetscInt off;

1129:       PetscDSGetComponentOffset(prob, f, &off);
1130:       if (isFE[f]) {
1131:         const PetscInt *cone;
1132:         PetscInt        comp, coneSize, faceLocL, faceLocR, ldof, rdof, d;

1134:         xL = xR = NULL;
1135:         DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1136:         DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1137:         DMPlexGetCone(dm, cells[0], &cone);
1138:         DMPlexGetConeSize(dm, cells[0], &coneSize);
1139:         for (faceLocL = 0; faceLocL < coneSize; ++faceLocL) if (cone[faceLocL] == face) break;
1140:         if (faceLocL == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[0]);
1141:         DMPlexGetCone(dm, cells[1], &cone);
1142:         DMPlexGetConeSize(dm, cells[1], &coneSize);
1143:         for (faceLocR = 0; faceLocR < coneSize; ++faceLocR) if (cone[faceLocR] == face) break;
1144:         if (faceLocR == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[1]);
1145:         /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
1146:         EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
1147:         if (rdof == ldof) {EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
1148:         else              {PetscSectionGetFieldComponents(section, f, &comp); for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
1149:         DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1150:         DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1151:       } else {
1152:         PetscFV  fv;
1153:         PetscInt numComp, c;

1155:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
1156:         PetscFVGetNumComponents(fv, &numComp);
1157:         DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
1158:         DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
1159:         if (dmGrad) {
1160:           PetscReal dxL[3], dxR[3];

1162:           DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
1163:           DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
1164:           DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
1165:           DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
1166:           for (c = 0; c < numComp; ++c) {
1167:             uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
1168:             uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
1169:           }
1170:         } else {
1171:           for (c = 0; c < numComp; ++c) {
1172:             uLl[iface*Nc+off+c] = xL[c];
1173:             uRl[iface*Nc+off+c] = xR[c];
1174:           }
1175:         }
1176:       }
1177:     }
1178:     ++iface;
1179:   }
1180:   VecRestoreArrayRead(locX, &x);
1181:   VecRestoreArrayRead(faceGeometry, &facegeom);
1182:   VecRestoreArrayRead(cellGeometry, &cellgeom);
1183:   if (locGrad) {
1184:     VecRestoreArrayRead(locGrad, &lgrad);
1185:   }
1186:   PetscFree(isFE);
1187:   return(0);
1188: }

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

1195:   Input Parameters:
1196: + dm     - The DM
1197: . fStart - The first face to include
1198: . fEnd   - The first face to exclude
1199: . locX   - A local vector with the solution fields
1200: . locX_t - A local vector with solution field time derivatives, or NULL
1201: . faceGeometry - A local vector with face geometry
1202: . cellGeometry - A local vector with cell geometry
1203: - locaGrad - A local vector with field gradients, or NULL

1205:   Output Parameters:
1206: + uL - The field values at the left side of the face
1207: - uR - The field values at the right side of the face

1209:   Level: developer

1211: .seealso: DMPlexGetFaceFields()
1212: @*/
1213: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1214: {

1218:   DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uL);
1219:   DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uR);
1220:   return(0);
1221: }

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

1228:   Input Parameters:
1229: + dm     - The DM
1230: . fStart - The first face to include
1231: . fEnd   - The first face to exclude
1232: . faceGeometry - A local vector with face geometry
1233: - cellGeometry - A local vector with cell geometry

1235:   Output Parameters:
1236: + fgeom - The extract the face centroid and normal
1237: - vol   - The cell volume

1239:   Level: developer

1241: .seealso: DMPlexGetCellFields()
1242: @*/
1243: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1244: {
1245:   DM                 dmFace, dmCell;
1246:   DMLabel            ghostLabel;
1247:   const PetscScalar *facegeom, *cellgeom;
1248:   PetscInt           dim, numFaces = fEnd - fStart, iface, face;
1249:   PetscErrorCode     ierr;

1257:   DMGetDimension(dm, &dim);
1258:   DMPlexGetLabel(dm, "ghost", &ghostLabel);
1259:   VecGetDM(faceGeometry, &dmFace);
1260:   VecGetArrayRead(faceGeometry, &facegeom);
1261:   VecGetDM(cellGeometry, &dmCell);
1262:   VecGetArrayRead(cellGeometry, &cellgeom);
1263:   PetscMalloc1(numFaces, fgeom);
1264:   DMGetWorkArray(dm, numFaces*2, PETSC_SCALAR, vol);
1265:   for (face = fStart, iface = 0; face < fEnd; ++face) {
1266:     const PetscInt        *cells;
1267:     const PetscFVFaceGeom *fg;
1268:     const PetscFVCellGeom *cgL, *cgR;
1269:     PetscFVFaceGeom       *fgeoml = *fgeom;
1270:     PetscReal             *voll   = *vol;
1271:     PetscInt               ghost, d;

1273:     DMLabelGetValue(ghostLabel, face, &ghost);
1274:     if (ghost >= 0) continue;
1275:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1276:     DMPlexGetSupport(dm, face, &cells);
1277:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1278:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1279:     for (d = 0; d < dim; ++d) {
1280:       fgeoml[iface].centroid[d] = fg->centroid[d];
1281:       fgeoml[iface].normal[d]   = fg->normal[d];
1282:     }
1283:     voll[iface*2+0] = cgL->volume;
1284:     voll[iface*2+1] = cgR->volume;
1285:     ++iface;
1286:   }
1287:   VecRestoreArrayRead(faceGeometry, &facegeom);
1288:   VecRestoreArrayRead(cellGeometry, &cellgeom);
1289:   return(0);
1290: }

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

1297:   Input Parameters:
1298: + dm     - The DM
1299: . fStart - The first face to include
1300: . fEnd   - The first face to exclude
1301: . faceGeometry - A local vector with face geometry
1302: - cellGeometry - A local vector with cell geometry

1304:   Output Parameters:
1305: + fgeom - The extract the face centroid and normal
1306: - vol   - The cell volume

1308:   Level: developer

1310: .seealso: DMPlexGetFaceFields()
1311: @*/
1312: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1313: {

1317:   PetscFree(*fgeom);
1318:   DMRestoreWorkArray(dm, 0, PETSC_REAL, vol);
1319:   return(0);
1320: }

1324: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
1325: {
1326:   DM                 dmFace, dmCell, dmGrad;
1327:   DMLabel            ghostLabel;
1328:   PetscDS            prob;
1329:   PetscFV            fvm;
1330:   PetscLimiter       lim;
1331:   const PetscScalar *facegeom, *cellgeom, *x;
1332:   PetscScalar       *gr;
1333:   PetscReal         *cellPhi;
1334:   PetscInt           dim, face, cell, totDim, cStart, cEnd, cEndInterior;
1335:   PetscErrorCode     ierr;

1338:   DMGetDimension(dm, &dim);
1339:   DMGetDS(dm, &prob);
1340:   PetscDSGetTotalDimension(prob, &totDim);
1341:   DMPlexGetLabel(dm, "ghost", &ghostLabel);
1342:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fvm);
1343:   PetscFVGetLimiter(fvm, &lim);
1344:   VecGetDM(faceGeometry, &dmFace);
1345:   VecGetArrayRead(faceGeometry, &facegeom);
1346:   VecGetDM(cellGeometry, &dmCell);
1347:   VecGetArrayRead(cellGeometry, &cellgeom);
1348:   VecGetArrayRead(locX, &x);
1349:   VecGetDM(grad, &dmGrad);
1350:   VecZeroEntries(grad);
1351:   VecGetArray(grad, &gr);
1352:   /* Reconstruct gradients */
1353:   for (face = fStart; face < fEnd; ++face) {
1354:     const PetscInt        *cells;
1355:     const PetscFVFaceGeom *fg;
1356:     const PetscScalar     *cx[2];
1357:     PetscScalar           *cgrad[2];
1358:     PetscBool              boundary;
1359:     PetscInt               ghost, c, pd, d;

1361:     DMLabelGetValue(ghostLabel, face, &ghost);
1362:     if (ghost >= 0) continue;
1363:     DMPlexIsBoundaryPoint(dm, face, &boundary);
1364:     if (boundary) continue;
1365:     DMPlexGetSupport(dm, face, &cells);
1366:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1367:     for (c = 0; c < 2; ++c) {
1368:       DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
1369:       DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]);
1370:     }
1371:     for (pd = 0; pd < totDim; ++pd) {
1372:       PetscScalar delta = cx[1][pd] - cx[0][pd];

1374:       for (d = 0; d < dim; ++d) {
1375:         if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta;
1376:         if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta;
1377:       }
1378:     }
1379:   }
1380:   /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
1381:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1382:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1383:   cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
1384:   DMGetWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1385:   for (cell = dmGrad && lim ? cStart : cEnd; cell < cEndInterior; ++cell) {
1386:     const PetscInt        *faces;
1387:     const PetscScalar     *cx;
1388:     const PetscFVCellGeom *cg;
1389:     PetscScalar           *cgrad;
1390:     PetscInt               coneSize, f, pd, d;

1392:     DMPlexGetConeSize(dm, cell, &coneSize);
1393:     DMPlexGetCone(dm, cell, &faces);
1394:     DMPlexPointLocalRead(dm, cell, x, &cx);
1395:     DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
1396:     DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
1397:     if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
1398:     /* Limiter will be minimum value over all neighbors */
1399:     for (d = 0; d < totDim; ++d) cellPhi[d] = PETSC_MAX_REAL;
1400:     for (f = 0; f < coneSize; ++f) {
1401:       const PetscScalar     *ncx;
1402:       const PetscFVCellGeom *ncg;
1403:       const PetscInt        *fcells;
1404:       PetscInt               face = faces[f], ncell, ghost;
1405:       PetscReal              v[3];
1406:       PetscBool              boundary;

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

1420:         PetscLimiterLimit(lim, flim, &phi);
1421:         cellPhi[d] = PetscMin(cellPhi[d], phi);
1422:       }
1423:     }
1424:     /* Apply limiter to gradient */
1425:     for (pd = 0; pd < totDim; ++pd)
1426:       /* Scalar limiter applied to each component separately */
1427:       for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd];
1428:   }
1429:   DMRestoreWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1430:   VecRestoreArrayRead(faceGeometry, &facegeom);
1431:   VecRestoreArrayRead(cellGeometry, &cellgeom);
1432:   VecRestoreArrayRead(locX, &x);
1433:   VecRestoreArray(grad, &gr);
1434:   return(0);
1435: }

1439: PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX_t, Vec locF, void *user)
1440: {
1441:   DM_Plex         *mesh = (DM_Plex *) dm->data;
1442:   PetscSection     section;
1443:   PetscDS          prob;
1444:   DMLabel          depth;
1445:   PetscFECellGeom *cgeom;
1446:   PetscScalar     *u = NULL, *u_t = NULL, *elemVec = NULL;
1447:   PetscInt         dim, Nf, f, totDimBd, numBd, bd;
1448:   PetscErrorCode   ierr;

1451:   DMGetDimension(dm, &dim);
1452:   DMGetDefaultSection(dm, &section);
1453:   DMGetDS(dm, &prob);
1454:   PetscDSGetNumFields(prob, &Nf);
1455:   PetscDSGetTotalBdDimension(prob, &totDimBd);
1456:   DMPlexGetDepthLabel(dm, &depth);
1457:   DMPlexGetNumBoundary(dm, &numBd);
1458:   for (bd = 0; bd < numBd; ++bd) {
1459:     const char     *bdLabel;
1460:     DMLabel         label;
1461:     IS              pointIS;
1462:     const PetscInt *points;
1463:     const PetscInt *values;
1464:     PetscInt        field, numValues, numPoints, p, dep, numFaces;
1465:     PetscBool       isEssential;

1467:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
1468:     if (isEssential) continue;
1469:     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1470:     DMPlexGetLabel(dm, bdLabel, &label);
1471:     DMLabelGetStratumSize(label, 1, &numPoints);
1472:     DMLabelGetStratumIS(label, 1, &pointIS);
1473:     ISGetIndices(pointIS, &points);
1474:     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1475:       DMLabelGetValue(depth, points[p], &dep);
1476:       if (dep == dim-1) ++numFaces;
1477:     }
1478:     PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd,&elemVec);
1479:     if (locX_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
1480:     for (p = 0, f = 0; p < numPoints; ++p) {
1481:       const PetscInt point = points[p];
1482:       PetscScalar   *x     = NULL;
1483:       PetscInt       i;

1485:       DMLabelGetValue(depth, points[p], &dep);
1486:       if (dep != dim-1) continue;
1487:       DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
1488:       DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
1489:       if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
1490:       DMPlexVecGetClosure(dm, section, locX, point, NULL, &x);
1491:       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1492:       DMPlexVecRestoreClosure(dm, section, locX, point, NULL, &x);
1493:       if (locX_t) {
1494:         DMPlexVecGetClosure(dm, section, locX_t, point, NULL, &x);
1495:         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1496:         DMPlexVecRestoreClosure(dm, section, locX_t, point, NULL, &x);
1497:       }
1498:       ++f;
1499:     }
1500:     for (f = 0; f < Nf; ++f) {
1501:       PetscFE         fe;
1502:       PetscQuadrature q;
1503:       PetscInt        numQuadPoints, Nb;
1504:       /* Conforming batches */
1505:       PetscInt        numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1506:       /* Remainder */
1507:       PetscInt        Nr, offset;

1509:       PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);
1510:       PetscFEGetQuadrature(fe, &q);
1511:       PetscFEGetDimension(fe, &Nb);
1512:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1513:       PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1514:       blockSize = Nb*numQuadPoints;
1515:       batchSize = numBlocks * blockSize;
1516:        PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1517:       numChunks = numFaces / (numBatches*batchSize);
1518:       Ne        = numChunks*numBatches*batchSize;
1519:       Nr        = numFaces % (numBatches*batchSize);
1520:       offset    = numFaces - Nr;
1521:       PetscFEIntegrateBdResidual(fe, prob, f, Ne, cgeom, u, u_t, NULL, NULL, elemVec);
1522:       PetscFEIntegrateBdResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);
1523:     }
1524:     for (p = 0, f = 0; p < numPoints; ++p) {
1525:       const PetscInt point = points[p];

1527:       DMLabelGetValue(depth, point, &dep);
1528:       if (dep != dim-1) continue;
1529:       if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);}
1530:       DMPlexVecSetClosure(dm, NULL, locF, point, &elemVec[f*totDimBd], ADD_VALUES);
1531:       ++f;
1532:     }
1533:     ISRestoreIndices(pointIS, &points);
1534:     ISDestroy(&pointIS);
1535:     PetscFree3(u,cgeom,elemVec);
1536:     if (locX_t) {PetscFree(u_t);}
1537:   }
1538:   return(0);
1539: }

1543: PetscErrorCode DMPlexComputeResidual_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
1544: {
1545:   DM_Plex          *mesh       = (DM_Plex *) dm->data;
1546:   const char       *name       = "Residual";
1547:   DM                dmAux      = NULL;
1548:   DM                dmGrad     = NULL;
1549:   DMLabel           ghostLabel = NULL;
1550:   PetscDS           prob       = NULL;
1551:   PetscDS           probAux    = NULL;
1552:   PetscSection      section    = NULL;
1553:   PetscBool         useFEM     = PETSC_FALSE;
1554:   PetscBool         useFVM     = PETSC_FALSE;
1555:   PetscBool         isImplicit = (locX_t || time == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
1556:   PetscFV           fvm        = NULL;
1557:   PetscFECellGeom  *cgeomFEM   = NULL;
1558:   PetscFVCellGeom  *cgeomFVM   = NULL;
1559:   PetscFVFaceGeom  *fgeomFVM   = NULL;
1560:   Vec               locA, cellGeometryFEM = NULL, cellGeometryFVM = NULL, faceGeometryFVM = NULL, grad, locGrad = NULL;
1561:   PetscScalar      *u, *u_t, *a, *uL, *uR;
1562:   PetscInt          Nf, f, totDim, totDimAux, numChunks, cellChunkSize, faceChunkSize, chunk, fStart, fEnd;
1563:   PetscErrorCode    ierr;

1566:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
1567:   /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
1568:   /* TODO The FVM geometry is over-manipulated. Make the precalc functions return exactly what we need */
1569:   /* FEM+FVM */
1570:   /* 1: Get sizes from dm and dmAux */
1571:   DMGetDefaultSection(dm, &section);
1572:   DMPlexGetLabel(dm, "ghost", &ghostLabel);
1573:   DMGetDS(dm, &prob);
1574:   PetscDSGetNumFields(prob, &Nf);
1575:   PetscDSGetTotalDimension(prob, &totDim);
1576:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1577:   if (locA) {
1578:     VecGetDM(locA, &dmAux);
1579:     DMGetDS(dmAux, &probAux);
1580:     PetscDSGetTotalDimension(probAux, &totDimAux);
1581:   }
1582:   /* 2: Get geometric data */
1583:   for (f = 0; f < Nf; ++f) {
1584:     PetscObject  obj;
1585:     PetscClassId id;
1586:     PetscBool    fimp;

1588:     PetscDSGetImplicit(prob, f, &fimp);
1589:     if (isImplicit != fimp) continue;
1590:     PetscDSGetDiscretization(prob, f, &obj);
1591:     PetscObjectGetClassId(obj, &id);
1592:     if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
1593:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1594:   }
1595:   if (useFEM) {
1596:     DMPlexSNESGetGeometryFEM(dm, &cellGeometryFEM);
1597:     VecGetArrayRead(cellGeometryFEM, (const PetscScalar **) &cgeomFEM);
1598:   }
1599:   if (useFVM) {
1600:     DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1601:     VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1602:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1603:     /* Reconstruct and limit cell gradients */
1604:     DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1605:     if (dmGrad) {
1606:       DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1607:       DMGetGlobalVector(dmGrad, &grad);
1608:       DMPlexReconstructGradients_Internal(dm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1609:       /* Communicate gradient values */
1610:       DMGetLocalVector(dmGrad, &locGrad);
1611:       DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1612:       DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1613:       DMRestoreGlobalVector(dmGrad, &grad);
1614:     }
1615:   }
1616:   /* Handle boundary values */
1617:   DMPlexInsertBoundaryValues(dm, locX, time, faceGeometryFVM, cellGeometryFVM, locGrad);
1618:   /* Loop over chunks */
1619:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1620:   numChunks     = 1;
1621:   cellChunkSize = (cEnd - cStart)/numChunks;
1622:   faceChunkSize = (fEnd - fStart)/numChunks;
1623:   for (chunk = 0; chunk < numChunks; ++chunk) {
1624:     PetscScalar     *elemVec, *fluxL, *fluxR;
1625:     PetscReal       *vol;
1626:     PetscFVFaceGeom *fgeom;
1627:     PetscInt         cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, cell;
1628:     PetscInt         fS = fStart+chunk*faceChunkSize, fE = PetscMin(fS+faceChunkSize, fEnd), numFaces = fE - fS, face;

1630:     /* Extract field coefficients */
1631:     if (useFEM) {
1632:       DMPlexGetCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1633:       DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1634:       PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
1635:     }
1636:     if (useFVM) {
1637:       DMPlexGetFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1638:       DMPlexGetFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1639:       DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1640:       DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1641:       PetscMemzero(fluxL, numFaces*totDim * sizeof(PetscScalar));
1642:       PetscMemzero(fluxR, numFaces*totDim * sizeof(PetscScalar));
1643:     }
1644:     /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
1645:     /* Loop over fields */
1646:     for (f = 0; f < Nf; ++f) {
1647:       PetscObject  obj;
1648:       PetscClassId id;
1649:       PetscBool    fimp;
1650:       PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1652:       PetscDSGetImplicit(prob, f, &fimp);
1653:       if (isImplicit != fimp) continue;
1654:       PetscDSGetDiscretization(prob, f, &obj);
1655:       PetscObjectGetClassId(obj, &id);
1656:       if (id == PETSCFE_CLASSID) {
1657:         PetscFE         fe = (PetscFE) obj;
1658:         PetscQuadrature q;
1659:         PetscInt        Nq, Nb;

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

1663:         PetscFEGetQuadrature(fe, &q);
1664:         PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1665:         PetscFEGetDimension(fe, &Nb);
1666:         blockSize = Nb*Nq;
1667:         batchSize = numBlocks * blockSize;
1668:          PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1669:         numChunks = numCells / (numBatches*batchSize);
1670:         Ne        = numChunks*numBatches*batchSize;
1671:         Nr        = numCells % (numBatches*batchSize);
1672:         offset    = numCells - Nr;
1673:         /* Integrate FE residual to get elemVec (need fields at quadrature points) */
1674:         /*   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) */
1675:         PetscFEIntegrateResidual(fe, prob, f, Ne, cgeomFEM, u, u_t, probAux, a, elemVec);
1676:         PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1677:       } else if (id == PETSCFV_CLASSID) {
1678:         PetscFV fv = (PetscFV) obj;

1680:         Ne = numFaces;
1681:         Nr = 0;
1682:         /* Riemann solve over faces (need fields at face centroids) */
1683:         /*   We need to evaluate FE fields at those coordinates */
1684:         PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
1685:       } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1686:     }
1687:     /* Loop over domain */
1688:     if (useFEM) {
1689:       /* Add elemVec to locX */
1690:       for (cell = cS; cell < cE; ++cell) {
1691:         if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cell*totDim]);}
1692:         DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cell*totDim], ADD_VALUES);
1693:       }
1694:     }
1695:     if (useFVM) {
1696:       PetscScalar *fa;
1697:       PetscInt     iface;

1699:       VecGetArray(locF, &fa);
1700:       for (f = 0; f < Nf; ++f) {
1701:         PetscFV      fv;
1702:         PetscObject  obj;
1703:         PetscClassId id;
1704:         PetscInt     foff, pdim;

1706:         PetscDSGetDiscretization(prob, f, &obj);
1707:         PetscDSGetFieldOffset(prob, f, &foff);
1708:         PetscObjectGetClassId(obj, &id);
1709:         if (id != PETSCFV_CLASSID) continue;
1710:         fv   = (PetscFV) obj;
1711:         PetscFVGetNumComponents(fv, &pdim);
1712:         /* Accumulate fluxes to cells */
1713:         for (face = fS, iface = 0; face < fE; ++face) {
1714:           const PetscInt *cells;
1715:           PetscScalar    *fL, *fR;
1716:           PetscInt        ghost, d;

1718:           DMLabelGetValue(ghostLabel, face, &ghost);
1719:           if (ghost >= 0) continue;
1720:           DMPlexGetSupport(dm, face, &cells);
1721:           DMPlexPointGlobalFieldRef(dm, cells[0], f, fa, &fL);
1722:           DMPlexPointGlobalFieldRef(dm, cells[1], f, fa, &fR);
1723:           for (d = 0; d < pdim; ++d) {
1724:             if (fL) fL[d] -= fluxL[iface*totDim+foff+d];
1725:             if (fR) fR[d] += fluxR[iface*totDim+foff+d];
1726:           }
1727:           ++iface;
1728:         }
1729:       }
1730:       VecRestoreArray(locF, &fa);
1731:     }
1732:     /* Handle time derivative */
1733:     if (locX_t) {
1734:       PetscScalar *x_t, *fa;

1736:       VecGetArray(locF, &fa);
1737:       VecGetArray(locX_t, &x_t);
1738:       for (f = 0; f < Nf; ++f) {
1739:         PetscFV      fv;
1740:         PetscObject  obj;
1741:         PetscClassId id;
1742:         PetscInt     pdim, d;

1744:         PetscDSGetDiscretization(prob, f, &obj);
1745:         PetscObjectGetClassId(obj, &id);
1746:         if (id != PETSCFV_CLASSID) continue;
1747:         fv   = (PetscFV) obj;
1748:         PetscFVGetNumComponents(fv, &pdim);
1749:         for (cell = cS; cell < cE; ++cell) {
1750:           PetscScalar *u_t, *r;

1752:           DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
1753:           DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
1754:           for (d = 0; d < pdim; ++d) r[d] += u_t[d];
1755:         }
1756:       }
1757:       VecRestoreArray(locX_t, &x_t);
1758:       VecRestoreArray(locF, &fa);
1759:     }
1760:     if (useFEM) {
1761:       DMPlexRestoreCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1762:       DMRestoreWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1763:     }
1764:     if (useFVM) {
1765:       DMPlexRestoreFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1766:       DMPlexRestoreFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1767:       DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1768:       DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1769:       if (dmGrad) {DMRestoreLocalVector(dmGrad, &locGrad);}
1770:     }
1771:   }

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

1775:   /* FEM */
1776:   /* 1: Get sizes from dm and dmAux */
1777:   /* 2: Get geometric data */
1778:   /* 3: Handle boundary values */
1779:   /* 4: Loop over domain */
1780:   /*   Extract coefficients */
1781:   /* Loop over fields */
1782:   /*   Set tiling for FE*/
1783:   /*   Integrate FE residual to get elemVec */
1784:   /*     Loop over subdomain */
1785:   /*       Loop over quad points */
1786:   /*         Transform coords to real space */
1787:   /*         Evaluate field and aux fields at point */
1788:   /*         Evaluate residual at point */
1789:   /*         Transform residual to real space */
1790:   /*       Add residual to elemVec */
1791:   /* Loop over domain */
1792:   /*   Add elemVec to locX */

1794:   /* FVM */
1795:   /* Get geometric data */
1796:   /* If using gradients */
1797:   /*   Compute gradient data */
1798:   /*   Loop over domain faces */
1799:   /*     Count computational faces */
1800:   /*     Reconstruct cell gradient */
1801:   /*   Loop over domain cells */
1802:   /*     Limit cell gradients */
1803:   /* Handle boundary values */
1804:   /* Loop over domain faces */
1805:   /*   Read out field, centroid, normal, volume for each side of face */
1806:   /* Riemann solve over faces */
1807:   /* Loop over domain faces */
1808:   /*   Accumulate fluxes to cells */
1809:   /* TODO Change printFEM to printDisc here */
1810:   if (mesh->printFEM) {DMPrintLocalVec(dm, name, mesh->printTol, locF);}
1811:   PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
1812:   return(0);
1813: }

1817: static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
1818: {
1819:   DM                dmCh, dmAux;
1820:   Vec               A, cellgeom;
1821:   PetscDS           prob, probCh, probAux = NULL;
1822:   PetscQuadrature   q;
1823:   PetscSection      section, sectionAux;
1824:   PetscFECellGeom  *cgeom;
1825:   PetscScalar      *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1826:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
1827:   PetscInt          totDim, totDimAux, diffCell = 0;
1828:   PetscErrorCode    ierr;

1831:   DMGetDimension(dm, &dim);
1832:   DMGetDefaultSection(dm, &section);
1833:   DMGetDS(dm, &prob);
1834:   PetscDSGetTotalDimension(prob, &totDim);
1835:   PetscSectionGetNumFields(section, &Nf);
1836:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1837:   numCells = cEnd - cStart;
1838:   PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);
1839:   DMGetDS(dmCh, &probCh);
1840:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1841:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1842:   if (dmAux) {
1843:     DMGetDefaultSection(dmAux, &sectionAux);
1844:     DMGetDS(dmAux, &probAux);
1845:     PetscDSGetTotalDimension(probAux, &totDimAux);
1846:   }
1847:   DMPlexInsertBoundaryValues(dm, X, 0.0, NULL, NULL, NULL);
1848:   VecSet(F, 0.0);
1849:   PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);
1850:   PetscMalloc1(numCells*totDim,&elemVecCh);
1851:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1852:   DMPlexSNESGetGeometryFEM(dm, &cellgeom);
1853:   VecGetArray(cellgeom, (PetscScalar **) &cgeom);
1854:   for (c = cStart; c < cEnd; ++c) {
1855:     PetscScalar *x = NULL, *x_t = NULL;
1856:     PetscInt     i;

1858:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
1859:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1860:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
1861:     if (X_t) {
1862:       DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
1863:       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1864:       DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
1865:     }
1866:     if (dmAux) {
1867:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1868:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1869:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1870:     }
1871:   }
1872:   for (f = 0; f < Nf; ++f) {
1873:     PetscFE  fe, feCh;
1874:     PetscInt numQuadPoints, Nb;
1875:     /* Conforming batches */
1876:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1877:     /* Remainder */
1878:     PetscInt Nr, offset;

1880:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1881:     PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);
1882:     PetscFEGetQuadrature(fe, &q);
1883:     PetscFEGetDimension(fe, &Nb);
1884:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1885:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1886:     blockSize = Nb*numQuadPoints;
1887:     batchSize = numBlocks * blockSize;
1888:      PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1889:     numChunks = numCells / (numBatches*batchSize);
1890:     Ne        = numChunks*numBatches*batchSize;
1891:     Nr        = numCells % (numBatches*batchSize);
1892:     offset    = numCells - Nr;
1893:     PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVec);
1894:     PetscFEIntegrateResidual(feCh, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVecCh);
1895:     PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1896:     PetscFEIntegrateResidual(feCh, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVecCh[offset*totDim]);
1897:   }
1898:   for (c = cStart; c < cEnd; ++c) {
1899:     PetscBool diff = PETSC_FALSE;
1900:     PetscInt  d;

1902:     for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1903:     if (diff) {
1904:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);
1905:       DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);
1906:       DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);
1907:       ++diffCell;
1908:     }
1909:     if (diffCell > 9) break;
1910:     DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);
1911:   }
1912:   VecRestoreArray(cellgeom, (PetscScalar **) &cgeom);
1913:   PetscFree3(u,u_t,elemVec);
1914:   PetscFree(elemVecCh);
1915:   if (dmAux) {PetscFree(a);}
1916:   return(0);
1917: }

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

1924:   Input Parameters:
1925: + dm - The mesh
1926: . X  - Local solution
1927: - user - The user context

1929:   Output Parameter:
1930: . F  - Local output vector

1932:   Level: developer

1934: .seealso: DMPlexComputeJacobianActionFEM()
1935: @*/
1936: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1937: {
1938:   PetscObject    check;
1939:   PetscInt       cStart, cEnd, cEndInterior;

1943:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1944:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1945:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1946:   /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
1947:   PetscObjectQuery((PetscObject) dm, "dmCh", &check);
1948:   if (check) {DMPlexComputeResidualFEM_Check_Internal(dm, X, NULL, F, user);}
1949:   else       {DMPlexComputeResidual_Internal(dm, cStart, cEnd, PETSC_MIN_REAL, X, NULL, F, user);}
1950:   return(0);
1951: }

1955: 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)
1956: {
1957:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1958:   const char       *name  = "Jacobian";
1959:   DM                dmAux;
1960:   DMLabel           depth;
1961:   Vec               A, cellgeom;
1962:   PetscDS           prob, probAux = NULL;
1963:   PetscQuadrature   quad;
1964:   PetscSection      section, globalSection, sectionAux;
1965:   PetscFECellGeom  *cgeom;
1966:   PetscScalar      *elemMat, *u, *u_t, *a = NULL;
1967:   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, c;
1968:   PetscInt          totDim, totDimBd, totDimAux, numBd, bd;
1969:   PetscBool         isShell;
1970:   PetscErrorCode    ierr;

1973:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
1974:   DMGetDimension(dm, &dim);
1975:   DMGetDefaultSection(dm, &section);
1976:   DMGetDefaultGlobalSection(dm, &globalSection);
1977:   DMGetDS(dm, &prob);
1978:   PetscDSGetTotalDimension(prob, &totDim);
1979:   PetscDSGetTotalBdDimension(prob, &totDimBd);
1980:   PetscSectionGetNumFields(section, &Nf);
1981:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1982:   numCells = cEnd - cStart;
1983:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1984:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1985:   if (dmAux) {
1986:     DMGetDefaultSection(dmAux, &sectionAux);
1987:     DMGetDS(dmAux, &probAux);
1988:     PetscDSGetTotalDimension(probAux, &totDimAux);
1989:   }
1990:   DMPlexInsertBoundaryValues(dm, X, 0.0, NULL, NULL, NULL);
1991:   MatZeroEntries(JacP);
1992:   PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat);
1993:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1994:   DMPlexSNESGetGeometryFEM(dm, &cellgeom);
1995:   VecGetArray(cellgeom, (PetscScalar **) &cgeom);
1996:   for (c = cStart; c < cEnd; ++c) {
1997:     PetscScalar *x = NULL,  *x_t = NULL;
1998:     PetscInt     i;

2000:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2001:     for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2002:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2003:     if (X_t) {
2004:       DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2005:       for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2006:       DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2007:     }
2008:     if (dmAux) {
2009:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
2010:       for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2011:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
2012:     }
2013:   }
2014:   PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));
2015:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
2016:     PetscFE  fe;
2017:     PetscInt numQuadPoints, Nb;
2018:     /* Conforming batches */
2019:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2020:     /* Remainder */
2021:     PetscInt Nr, offset;

2023:     PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2024:     PetscFEGetQuadrature(fe, &quad);
2025:     PetscFEGetDimension(fe, &Nb);
2026:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2027:     PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2028:     blockSize = Nb*numQuadPoints;
2029:     batchSize = numBlocks * blockSize;
2030:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2031:     numChunks = numCells / (numBatches*batchSize);
2032:     Ne        = numChunks*numBatches*batchSize;
2033:     Nr        = numCells % (numBatches*batchSize);
2034:     offset    = numCells - Nr;
2035:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2036:       PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMat);
2037:       PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMat[offset*totDim*totDim]);
2038:     }
2039:   }
2040:   for (c = cStart; c < cEnd; ++c) {
2041:     if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2042:     DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2043:   }
2044:   VecGetArray(cellgeom, (PetscScalar **) &cgeom);
2045:   PetscFree3(u,u_t,elemMat);
2046:   if (dmAux) {PetscFree(a);}
2047:   DMPlexGetDepthLabel(dm, &depth);
2048:   DMPlexGetNumBoundary(dm, &numBd);
2049:   DMPlexGetDepthLabel(dm, &depth);
2050:   DMPlexGetNumBoundary(dm, &numBd);
2051:   for (bd = 0; bd < numBd; ++bd) {
2052:     const char     *bdLabel;
2053:     DMLabel         label;
2054:     IS              pointIS;
2055:     const PetscInt *points;
2056:     const PetscInt *values;
2057:     PetscInt        field, numValues, numPoints, p, dep, numFaces;
2058:     PetscBool       isEssential;

2060:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
2061:     if (isEssential) continue;
2062:     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
2063:     DMPlexGetLabel(dm, bdLabel, &label);
2064:     DMLabelGetStratumSize(label, 1, &numPoints);
2065:     DMLabelGetStratumIS(label, 1, &pointIS);
2066:     ISGetIndices(pointIS, &points);
2067:     for (p = 0, numFaces = 0; p < numPoints; ++p) {
2068:       DMLabelGetValue(depth, points[p], &dep);
2069:       if (dep == dim-1) ++numFaces;
2070:     }
2071:     PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd*totDimBd,&elemMat);
2072:     if (X_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
2073:     for (p = 0, f = 0; p < numPoints; ++p) {
2074:       const PetscInt point = points[p];
2075:       PetscScalar   *x     = NULL;
2076:       PetscInt       i;

2078:       DMLabelGetValue(depth, points[p], &dep);
2079:       if (dep != dim-1) continue;
2080:       DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
2081:       DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
2082:       if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
2083:       DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
2084:       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
2085:       DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
2086:       if (X_t) {
2087:         DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);
2088:         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
2089:         DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);
2090:       }
2091:       ++f;
2092:     }
2093:     PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));
2094:     for (fieldI = 0; fieldI < Nf; ++fieldI) {
2095:       PetscFE  fe;
2096:       PetscInt numQuadPoints, Nb;
2097:       /* Conforming batches */
2098:       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2099:       /* Remainder */
2100:       PetscInt Nr, offset;

2102:       PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
2103:       PetscFEGetQuadrature(fe, &quad);
2104:       PetscFEGetDimension(fe, &Nb);
2105:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2106:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2107:       blockSize = Nb*numQuadPoints;
2108:       batchSize = numBlocks * blockSize;
2109:        PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2110:       numChunks = numFaces / (numBatches*batchSize);
2111:       Ne        = numChunks*numBatches*batchSize;
2112:       Nr        = numFaces % (numBatches*batchSize);
2113:       offset    = numFaces - Nr;
2114:       for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2115:         PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, NULL, NULL, elemMat);
2116:         PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);
2117:       }
2118:     }
2119:     for (p = 0, f = 0; p < numPoints; ++p) {
2120:       const PetscInt point = points[p];

2122:       DMLabelGetValue(depth, point, &dep);
2123:       if (dep != dim-1) continue;
2124:       if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);}
2125:       DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);
2126:       ++f;
2127:     }
2128:     ISRestoreIndices(pointIS, &points);
2129:     ISDestroy(&pointIS);
2130:     PetscFree3(u,cgeom,elemMat);
2131:     if (X_t) {PetscFree(u_t);}
2132:   }
2133:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
2134:   MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
2135:   if (mesh->printFEM) {
2136:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
2137:     MatChop(JacP, 1.0e-10);
2138:     MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
2139:   }
2140:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2141:   PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
2142:   if (isShell) {
2143:     JacActionCtx *jctx;

2145:     MatShellGetContext(Jac, &jctx);
2146:     VecCopy(X, jctx->u);
2147:   }
2148:   return(0);
2149: }

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

2156:   Input Parameters:
2157: + dm - The mesh
2158: . X  - Local input vector
2159: - user - The user context

2161:   Output Parameter:
2162: . Jac  - Jacobian matrix

2164:   Note:
2165:   The first member of the user context must be an FEMContext.

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

2170:   Level: developer

2172: .seealso: FormFunctionLocal()
2173: @*/
2174: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
2175: {
2176:   PetscInt       cStart, cEnd, cEndInterior;

2180:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2181:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2182:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2183:   DMPlexComputeJacobian_Internal(dm, cStart, cEnd, 0.0, 0.0, X, NULL, Jac, JacP, user);
2184:   return(0);
2185: }

2189: PetscErrorCode DMSNESCheckFromOptions_Internal(SNES snes, DM dm, Vec u, Vec sol, PetscErrorCode (**exactFuncs)(PetscInt, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx), void **ctxs)
2190: {
2191:   Mat            J, M;
2192:   Vec            r, b;
2193:   MatNullSpace   nullSpace;
2194:   PetscReal     *error, res = 0.0;
2195:   PetscInt       numFields;

2199:   VecDuplicate(u, &r);
2200:   DMCreateMatrix(dm, &J);
2201:   M    = J;
2202:   /* TODO Null space for J */
2203:   /* Check discretization error */
2204:   DMGetNumFields(dm, &numFields);
2205:   PetscMalloc1(PetscMax(1, numFields), &error);
2206:   if (numFields > 1) {
2207:     PetscInt f;

2209:     DMPlexComputeL2FieldDiff(dm, exactFuncs, ctxs, u, error);
2210:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: [");
2211:     for (f = 0; f < numFields; ++f) {
2212:       if (f) {PetscPrintf(PETSC_COMM_WORLD, ", ");}
2213:       if (error[f] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "%g", error[f]);}
2214:       else                     {PetscPrintf(PETSC_COMM_WORLD, "< 1.0e-11");}
2215:     }
2216:     PetscPrintf(PETSC_COMM_WORLD, "]\n");
2217:   } else {
2218:     DMPlexComputeL2Diff(dm, exactFuncs, ctxs, u, &error[0]);
2219:     if (error[0] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error[0]);}
2220:     else                     {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
2221:   }
2222:   PetscFree(error);
2223:   /* Check residual */
2224:   SNESComputeFunction(snes, u, r);
2225:   VecNorm(r, NORM_2, &res);
2226:   PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
2227:   VecChop(r, 1.0e-10);
2228:   PetscObjectSetName((PetscObject) r, "Initial Residual");
2229:   PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
2230:   VecViewFromOptions(r, NULL, "-vec_view");
2231:   /* Check Jacobian */
2232:   SNESComputeJacobian(snes, u, M, M);
2233:   MatGetNullSpace(J, &nullSpace);
2234:   if (nullSpace) {
2235:     PetscBool isNull;
2236:     MatNullSpaceTest(nullSpace, J, &isNull);
2237:     if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
2238:   }
2239:   VecDuplicate(u, &b);
2240:   VecSet(r, 0.0);
2241:   SNESComputeFunction(snes, r, b);
2242:   MatMult(M, u, r);
2243:   VecAXPY(r, 1.0, b);
2244:   VecDestroy(&b);
2245:   VecNorm(r, NORM_2, &res);
2246:   PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
2247:   VecChop(r, 1.0e-10);
2248:   PetscObjectSetName((PetscObject) r, "Au - b = Au + F(0)");
2249:   PetscObjectSetOptionsPrefix((PetscObject)r,"linear_res_");
2250:   VecViewFromOptions(r, NULL, "-vec_view");
2251:   VecDestroy(&r);
2252:   MatNullSpaceDestroy(&nullSpace);
2253:   MatDestroy(&J);
2254:   return(0);
2255: }

2259: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
2260: {
2261:   DM             dm;
2262:   Vec            sol;
2263:   PetscBool      check;

2267:   PetscOptionsHasName(snes->hdr.prefix, "-dmsnes_check", &check);
2268:   if (!check) return(0);
2269:   SNESGetDM(snes, &dm);
2270:   VecDuplicate(u, &sol);
2271:   SNESSetSolution(snes, sol);
2272:   DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
2273:   VecDestroy(&sol);
2274:   return(0);
2275: }