xref: /petsc/src/snes/utils/dmplexsnes.c (revision c49ccbb32729e871fff84de9e1f3d01f6ca029ee) 
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 
6 /************************** Interpolation *******************************/
7 
8 #undef __FUNCT__
9 #define __FUNCT__ "DMInterpolationCreate"
10 PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
11 {
12   PetscErrorCode ierr;
13 
14   PetscFunctionBegin;
15   PetscValidPointer(ctx, 2);
16   ierr = PetscMalloc(sizeof(struct _DMInterpolationInfo), ctx);CHKERRQ(ierr);
17 
18   (*ctx)->comm   = comm;
19   (*ctx)->dim    = -1;
20   (*ctx)->nInput = 0;
21   (*ctx)->points = NULL;
22   (*ctx)->cells  = NULL;
23   (*ctx)->n      = -1;
24   (*ctx)->coords = NULL;
25   PetscFunctionReturn(0);
26 }
27 
28 #undef __FUNCT__
29 #define __FUNCT__ "DMInterpolationSetDim"
30 PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
31 {
32   PetscFunctionBegin;
33   if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);
34   ctx->dim = dim;
35   PetscFunctionReturn(0);
36 }
37 
38 #undef __FUNCT__
39 #define __FUNCT__ "DMInterpolationGetDim"
40 PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
41 {
42   PetscFunctionBegin;
43   PetscValidIntPointer(dim, 2);
44   *dim = ctx->dim;
45   PetscFunctionReturn(0);
46 }
47 
48 #undef __FUNCT__
49 #define __FUNCT__ "DMInterpolationSetDof"
50 PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
51 {
52   PetscFunctionBegin;
53   if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);
54   ctx->dof = dof;
55   PetscFunctionReturn(0);
56 }
57 
58 #undef __FUNCT__
59 #define __FUNCT__ "DMInterpolationGetDof"
60 PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
61 {
62   PetscFunctionBegin;
63   PetscValidIntPointer(dof, 2);
64   *dof = ctx->dof;
65   PetscFunctionReturn(0);
66 }
67 
68 #undef __FUNCT__
69 #define __FUNCT__ "DMInterpolationAddPoints"
70 PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
71 {
72   PetscErrorCode ierr;
73 
74   PetscFunctionBegin;
75   if (ctx->dim < 0) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
76   if (ctx->points)  SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
77   ctx->nInput = n;
78 
79   ierr = PetscMalloc1(n*ctx->dim, &ctx->points);CHKERRQ(ierr);
80   ierr = PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));CHKERRQ(ierr);
81   PetscFunctionReturn(0);
82 }
83 
84 #undef __FUNCT__
85 #define __FUNCT__ "DMInterpolationSetUp"
86 PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
87 {
88   MPI_Comm       comm = ctx->comm;
89   PetscScalar    *a;
90   PetscInt       p, q, i;
91   PetscMPIInt    rank, size;
92   PetscErrorCode ierr;
93   Vec            pointVec;
94   IS             cellIS;
95   PetscLayout    layout;
96   PetscReal      *globalPoints;
97   PetscScalar    *globalPointsScalar;
98   const PetscInt *ranges;
99   PetscMPIInt    *counts, *displs;
100   const PetscInt *foundCells;
101   PetscMPIInt    *foundProcs, *globalProcs;
102   PetscInt       n, N;
103 
104   PetscFunctionBegin;
105   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
106   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
107   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
108   if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
109   /* Locate points */
110   n = ctx->nInput;
111   if (!redundantPoints) {
112     ierr = PetscLayoutCreate(comm, &layout);CHKERRQ(ierr);
113     ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr);
114     ierr = PetscLayoutSetLocalSize(layout, n);CHKERRQ(ierr);
115     ierr = PetscLayoutSetUp(layout);CHKERRQ(ierr);
116     ierr = PetscLayoutGetSize(layout, &N);CHKERRQ(ierr);
117     /* Communicate all points to all processes */
118     ierr = PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);CHKERRQ(ierr);
119     ierr = PetscLayoutGetRanges(layout, &ranges);CHKERRQ(ierr);
120     for (p = 0; p < size; ++p) {
121       counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
122       displs[p] = ranges[p]*ctx->dim;
123     }
124     ierr = MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);CHKERRQ(ierr);
125   } else {
126     N = n;
127     globalPoints = ctx->points;
128     counts = displs = NULL;
129     layout = NULL;
130   }
131 #if 0
132   ierr = PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);CHKERRQ(ierr);
133   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
134 #else
135 #if defined(PETSC_USE_COMPLEX)
136   ierr = PetscMalloc1(N,&globalPointsScalar);CHKERRQ(ierr);
137   for (i=0; i<N; i++) globalPointsScalar[i] = globalPoints[i];
138 #else
139   globalPointsScalar = globalPoints;
140 #endif
141   ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);CHKERRQ(ierr);
142   ierr = PetscMalloc2(N,&foundProcs,N,&globalProcs);CHKERRQ(ierr);
143   ierr = DMLocatePoints(dm, pointVec, &cellIS);CHKERRQ(ierr);
144   ierr = ISGetIndices(cellIS, &foundCells);CHKERRQ(ierr);
145 #endif
146   for (p = 0; p < N; ++p) {
147     if (foundCells[p] >= 0) foundProcs[p] = rank;
148     else foundProcs[p] = size;
149   }
150   /* Let the lowest rank process own each point */
151   ierr   = MPI_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);CHKERRQ(ierr);
152   ctx->n = 0;
153   for (p = 0; p < N; ++p) {
154     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);
155     else if (globalProcs[p] == rank) ctx->n++;
156   }
157   /* Create coordinates vector and array of owned cells */
158   ierr = PetscMalloc1(ctx->n, &ctx->cells);CHKERRQ(ierr);
159   ierr = VecCreate(comm, &ctx->coords);CHKERRQ(ierr);
160   ierr = VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);CHKERRQ(ierr);
161   ierr = VecSetBlockSize(ctx->coords, ctx->dim);CHKERRQ(ierr);
162   ierr = VecSetType(ctx->coords,VECSTANDARD);CHKERRQ(ierr);
163   ierr = VecGetArray(ctx->coords, &a);CHKERRQ(ierr);
164   for (p = 0, q = 0, i = 0; p < N; ++p) {
165     if (globalProcs[p] == rank) {
166       PetscInt d;
167 
168       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
169       ctx->cells[q++] = foundCells[p];
170     }
171   }
172   ierr = VecRestoreArray(ctx->coords, &a);CHKERRQ(ierr);
173 #if 0
174   ierr = PetscFree3(foundCells,foundProcs,globalProcs);CHKERRQ(ierr);
175 #else
176   ierr = PetscFree2(foundProcs,globalProcs);CHKERRQ(ierr);
177   ierr = ISRestoreIndices(cellIS, &foundCells);CHKERRQ(ierr);
178   ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
179   ierr = VecDestroy(&pointVec);CHKERRQ(ierr);
180 #endif
181   if ((void*)globalPointsScalar != (void*)globalPoints) {ierr = PetscFree(globalPointsScalar);CHKERRQ(ierr);}
182   if (!redundantPoints) {ierr = PetscFree3(globalPoints,counts,displs);CHKERRQ(ierr);}
183   ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr);
184   PetscFunctionReturn(0);
185 }
186 
187 #undef __FUNCT__
188 #define __FUNCT__ "DMInterpolationGetCoordinates"
189 PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
190 {
191   PetscFunctionBegin;
192   PetscValidPointer(coordinates, 2);
193   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
194   *coordinates = ctx->coords;
195   PetscFunctionReturn(0);
196 }
197 
198 #undef __FUNCT__
199 #define __FUNCT__ "DMInterpolationGetVector"
200 PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
201 {
202   PetscErrorCode ierr;
203 
204   PetscFunctionBegin;
205   PetscValidPointer(v, 2);
206   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
207   ierr = VecCreate(ctx->comm, v);CHKERRQ(ierr);
208   ierr = VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);CHKERRQ(ierr);
209   ierr = VecSetBlockSize(*v, ctx->dof);CHKERRQ(ierr);
210   ierr = VecSetType(*v,VECSTANDARD);CHKERRQ(ierr);
211   PetscFunctionReturn(0);
212 }
213 
214 #undef __FUNCT__
215 #define __FUNCT__ "DMInterpolationRestoreVector"
216 PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
217 {
218   PetscErrorCode ierr;
219 
220   PetscFunctionBegin;
221   PetscValidPointer(v, 2);
222   if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
223   ierr = VecDestroy(v);CHKERRQ(ierr);
224   PetscFunctionReturn(0);
225 }
226 
227 #undef __FUNCT__
228 #define __FUNCT__ "DMInterpolate_Triangle_Private"
229 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
230 {
231   PetscReal      *v0, *J, *invJ, detJ;
232   PetscScalar    *a, *coords;
233   PetscInt       p;
234   PetscErrorCode ierr;
235 
236   PetscFunctionBegin;
237   ierr = PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);CHKERRQ(ierr);
238   ierr = VecGetArray(ctx->coords, &coords);CHKERRQ(ierr);
239   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
240   for (p = 0; p < ctx->n; ++p) {
241     PetscInt     c = ctx->cells[p];
242     PetscScalar *x = NULL;
243     PetscReal    xi[4];
244     PetscInt     d, f, comp;
245 
246     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
247     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
248     ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
249     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
250 
251     for (d = 0; d < ctx->dim; ++d) {
252       xi[d] = 0.0;
253       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
254       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];
255     }
256     ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
257   }
258   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
259   ierr = VecRestoreArray(ctx->coords, &coords);CHKERRQ(ierr);
260   ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr);
261   PetscFunctionReturn(0);
262 }
263 
264 #undef __FUNCT__
265 #define __FUNCT__ "DMInterpolate_Tetrahedron_Private"
266 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
267 {
268   PetscReal      *v0, *J, *invJ, detJ;
269   PetscScalar    *a, *coords;
270   PetscInt       p;
271   PetscErrorCode ierr;
272 
273   PetscFunctionBegin;
274   ierr = PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);CHKERRQ(ierr);
275   ierr = VecGetArray(ctx->coords, &coords);CHKERRQ(ierr);
276   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
277   for (p = 0; p < ctx->n; ++p) {
278     PetscInt       c = ctx->cells[p];
279     const PetscInt order[3] = {2, 1, 3};
280     PetscScalar   *x = NULL;
281     PetscReal      xi[4];
282     PetscInt       d, f, comp;
283 
284     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
285     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
286     ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
287     for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
288 
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[order[d]*ctx->dof+comp] - x[0*ctx->dof+comp])*xi[d];
293     }
294     ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
295   }
296   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
297   ierr = VecRestoreArray(ctx->coords, &coords);CHKERRQ(ierr);
298   ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr);
299   PetscFunctionReturn(0);
300 }
301 
302 #undef __FUNCT__
303 #define __FUNCT__ "QuadMap_Private"
304 PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
305 {
306   const PetscScalar *vertices = (const PetscScalar*) ctx;
307   const PetscScalar x0        = vertices[0];
308   const PetscScalar y0        = vertices[1];
309   const PetscScalar x1        = vertices[2];
310   const PetscScalar y1        = vertices[3];
311   const PetscScalar x2        = vertices[4];
312   const PetscScalar y2        = vertices[5];
313   const PetscScalar x3        = vertices[6];
314   const PetscScalar y3        = vertices[7];
315   const PetscScalar f_1       = x1 - x0;
316   const PetscScalar g_1       = y1 - y0;
317   const PetscScalar f_3       = x3 - x0;
318   const PetscScalar g_3       = y3 - y0;
319   const PetscScalar f_01      = x2 - x1 - x3 + x0;
320   const PetscScalar g_01      = y2 - y1 - y3 + y0;
321   PetscScalar       *ref, *real;
322   PetscErrorCode    ierr;
323 
324   PetscFunctionBegin;
325   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
326   ierr = VecGetArray(Xreal, &real);CHKERRQ(ierr);
327   {
328     const PetscScalar p0 = ref[0];
329     const PetscScalar p1 = ref[1];
330 
331     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
332     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
333   }
334   ierr = PetscLogFlops(28);CHKERRQ(ierr);
335   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
336   ierr = VecRestoreArray(Xreal, &real);CHKERRQ(ierr);
337   PetscFunctionReturn(0);
338 }
339 
340 #include <petsc-private/dmimpl.h>
341 #undef __FUNCT__
342 #define __FUNCT__ "QuadJacobian_Private"
343 PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, 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_01      = x2 - x1 - x3 + x0;
355   const PetscScalar g_01      = y2 - y1 - y3 + y0;
356   PetscScalar       *ref;
357   PetscErrorCode    ierr;
358 
359   PetscFunctionBegin;
360   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
361   {
362     const PetscScalar x       = ref[0];
363     const PetscScalar y       = ref[1];
364     const PetscInt    rows[2] = {0, 1};
365     PetscScalar       values[4];
366 
367     values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
368     values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
369     ierr      = MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);CHKERRQ(ierr);
370   }
371   ierr = PetscLogFlops(30);CHKERRQ(ierr);
372   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
373   ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
374   ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
375   PetscFunctionReturn(0);
376 }
377 
378 #undef __FUNCT__
379 #define __FUNCT__ "DMInterpolate_Quad_Private"
380 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
381 {
382   DM             dmCoord;
383   SNES           snes;
384   KSP            ksp;
385   PC             pc;
386   Vec            coordsLocal, r, ref, real;
387   Mat            J;
388   PetscScalar    *a, *coords;
389   PetscInt       p;
390   PetscErrorCode ierr;
391 
392   PetscFunctionBegin;
393   ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
394   ierr = DMGetCoordinateDM(dm, &dmCoord);CHKERRQ(ierr);
395   ierr = SNESCreate(PETSC_COMM_SELF, &snes);CHKERRQ(ierr);
396   ierr = SNESSetOptionsPrefix(snes, "quad_interp_");CHKERRQ(ierr);
397   ierr = VecCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
398   ierr = VecSetSizes(r, 2, 2);CHKERRQ(ierr);
399   ierr = VecSetType(r,dm->vectype);CHKERRQ(ierr);
400   ierr = VecDuplicate(r, &ref);CHKERRQ(ierr);
401   ierr = VecDuplicate(r, &real);CHKERRQ(ierr);
402   ierr = MatCreate(PETSC_COMM_SELF, &J);CHKERRQ(ierr);
403   ierr = MatSetSizes(J, 2, 2, 2, 2);CHKERRQ(ierr);
404   ierr = MatSetType(J, MATSEQDENSE);CHKERRQ(ierr);
405   ierr = MatSetUp(J);CHKERRQ(ierr);
406   ierr = SNESSetFunction(snes, r, QuadMap_Private, NULL);CHKERRQ(ierr);
407   ierr = SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);CHKERRQ(ierr);
408   ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
409   ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
410   ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
411   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
412 
413   ierr = VecGetArray(ctx->coords, &coords);CHKERRQ(ierr);
414   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
415   for (p = 0; p < ctx->n; ++p) {
416     PetscScalar *x = NULL, *vertices = NULL;
417     PetscScalar *xi;
418     PetscReal    xir[2];
419     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;
420 
421     /* Can make this do all points at once */
422     ierr = DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
423     if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);
424     ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
425     if (4*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 4*ctx->dof);
426     ierr   = SNESSetFunction(snes, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
427     ierr   = SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
428     ierr   = VecGetArray(real, &xi);CHKERRQ(ierr);
429     xi[0]  = coords[p*ctx->dim+0];
430     xi[1]  = coords[p*ctx->dim+1];
431     ierr   = VecRestoreArray(real, &xi);CHKERRQ(ierr);
432     ierr   = SNESSolve(snes, real, ref);CHKERRQ(ierr);
433     ierr   = VecGetArray(ref, &xi);CHKERRQ(ierr);
434     xir[0] = PetscRealPart(xi[0]);
435     xir[1] = PetscRealPart(xi[1]);
436     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];
437 
438     ierr = VecRestoreArray(ref, &xi);CHKERRQ(ierr);
439     ierr = DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
440     ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
441   }
442   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
443   ierr = VecRestoreArray(ctx->coords, &coords);CHKERRQ(ierr);
444 
445   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
446   ierr = VecDestroy(&r);CHKERRQ(ierr);
447   ierr = VecDestroy(&ref);CHKERRQ(ierr);
448   ierr = VecDestroy(&real);CHKERRQ(ierr);
449   ierr = MatDestroy(&J);CHKERRQ(ierr);
450   PetscFunctionReturn(0);
451 }
452 
453 #undef __FUNCT__
454 #define __FUNCT__ "HexMap_Private"
455 PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
456 {
457   const PetscScalar *vertices = (const PetscScalar*) ctx;
458   const PetscScalar x0        = vertices[0];
459   const PetscScalar y0        = vertices[1];
460   const PetscScalar z0        = vertices[2];
461   const PetscScalar x1        = vertices[9];
462   const PetscScalar y1        = vertices[10];
463   const PetscScalar z1        = vertices[11];
464   const PetscScalar x2        = vertices[6];
465   const PetscScalar y2        = vertices[7];
466   const PetscScalar z2        = vertices[8];
467   const PetscScalar x3        = vertices[3];
468   const PetscScalar y3        = vertices[4];
469   const PetscScalar z3        = vertices[5];
470   const PetscScalar x4        = vertices[12];
471   const PetscScalar y4        = vertices[13];
472   const PetscScalar z4        = vertices[14];
473   const PetscScalar x5        = vertices[15];
474   const PetscScalar y5        = vertices[16];
475   const PetscScalar z5        = vertices[17];
476   const PetscScalar x6        = vertices[18];
477   const PetscScalar y6        = vertices[19];
478   const PetscScalar z6        = vertices[20];
479   const PetscScalar x7        = vertices[21];
480   const PetscScalar y7        = vertices[22];
481   const PetscScalar z7        = vertices[23];
482   const PetscScalar f_1       = x1 - x0;
483   const PetscScalar g_1       = y1 - y0;
484   const PetscScalar h_1       = z1 - z0;
485   const PetscScalar f_3       = x3 - x0;
486   const PetscScalar g_3       = y3 - y0;
487   const PetscScalar h_3       = z3 - z0;
488   const PetscScalar f_4       = x4 - x0;
489   const PetscScalar g_4       = y4 - y0;
490   const PetscScalar h_4       = z4 - z0;
491   const PetscScalar f_01      = x2 - x1 - x3 + x0;
492   const PetscScalar g_01      = y2 - y1 - y3 + y0;
493   const PetscScalar h_01      = z2 - z1 - z3 + z0;
494   const PetscScalar f_12      = x7 - x3 - x4 + x0;
495   const PetscScalar g_12      = y7 - y3 - y4 + y0;
496   const PetscScalar h_12      = z7 - z3 - z4 + z0;
497   const PetscScalar f_02      = x5 - x1 - x4 + x0;
498   const PetscScalar g_02      = y5 - y1 - y4 + y0;
499   const PetscScalar h_02      = z5 - z1 - z4 + z0;
500   const PetscScalar f_012     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
501   const PetscScalar g_012     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
502   const PetscScalar h_012     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
503   PetscScalar       *ref, *real;
504   PetscErrorCode    ierr;
505 
506   PetscFunctionBegin;
507   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
508   ierr = VecGetArray(Xreal, &real);CHKERRQ(ierr);
509   {
510     const PetscScalar p0 = ref[0];
511     const PetscScalar p1 = ref[1];
512     const PetscScalar p2 = ref[2];
513 
514     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;
515     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;
516     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;
517   }
518   ierr = PetscLogFlops(114);CHKERRQ(ierr);
519   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
520   ierr = VecRestoreArray(Xreal, &real);CHKERRQ(ierr);
521   PetscFunctionReturn(0);
522 }
523 
524 #undef __FUNCT__
525 #define __FUNCT__ "HexJacobian_Private"
526 PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
527 {
528   const PetscScalar *vertices = (const PetscScalar*) ctx;
529   const PetscScalar x0        = vertices[0];
530   const PetscScalar y0        = vertices[1];
531   const PetscScalar z0        = vertices[2];
532   const PetscScalar x1        = vertices[9];
533   const PetscScalar y1        = vertices[10];
534   const PetscScalar z1        = vertices[11];
535   const PetscScalar x2        = vertices[6];
536   const PetscScalar y2        = vertices[7];
537   const PetscScalar z2        = vertices[8];
538   const PetscScalar x3        = vertices[3];
539   const PetscScalar y3        = vertices[4];
540   const PetscScalar z3        = vertices[5];
541   const PetscScalar x4        = vertices[12];
542   const PetscScalar y4        = vertices[13];
543   const PetscScalar z4        = vertices[14];
544   const PetscScalar x5        = vertices[15];
545   const PetscScalar y5        = vertices[16];
546   const PetscScalar z5        = vertices[17];
547   const PetscScalar x6        = vertices[18];
548   const PetscScalar y6        = vertices[19];
549   const PetscScalar z6        = vertices[20];
550   const PetscScalar x7        = vertices[21];
551   const PetscScalar y7        = vertices[22];
552   const PetscScalar z7        = vertices[23];
553   const PetscScalar f_xy      = x2 - x1 - x3 + x0;
554   const PetscScalar g_xy      = y2 - y1 - y3 + y0;
555   const PetscScalar h_xy      = z2 - z1 - z3 + z0;
556   const PetscScalar f_yz      = x7 - x3 - x4 + x0;
557   const PetscScalar g_yz      = y7 - y3 - y4 + y0;
558   const PetscScalar h_yz      = z7 - z3 - z4 + z0;
559   const PetscScalar f_xz      = x5 - x1 - x4 + x0;
560   const PetscScalar g_xz      = y5 - y1 - y4 + y0;
561   const PetscScalar h_xz      = z5 - z1 - z4 + z0;
562   const PetscScalar f_xyz     = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
563   const PetscScalar g_xyz     = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
564   const PetscScalar h_xyz     = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
565   PetscScalar       *ref;
566   PetscErrorCode    ierr;
567 
568   PetscFunctionBegin;
569   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
570   {
571     const PetscScalar x       = ref[0];
572     const PetscScalar y       = ref[1];
573     const PetscScalar z       = ref[2];
574     const PetscInt    rows[3] = {0, 1, 2};
575     PetscScalar       values[9];
576 
577     values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
578     values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
579     values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
580     values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
581     values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
582     values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
583     values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
584     values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
585     values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
586 
587     ierr = MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);CHKERRQ(ierr);
588   }
589   ierr = PetscLogFlops(152);CHKERRQ(ierr);
590   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
591   ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
592   ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
593   PetscFunctionReturn(0);
594 }
595 
596 #undef __FUNCT__
597 #define __FUNCT__ "DMInterpolate_Hex_Private"
598 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
599 {
600   DM             dmCoord;
601   SNES           snes;
602   KSP            ksp;
603   PC             pc;
604   Vec            coordsLocal, r, ref, real;
605   Mat            J;
606   PetscScalar    *a, *coords;
607   PetscInt       p;
608   PetscErrorCode ierr;
609 
610   PetscFunctionBegin;
611   ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
612   ierr = DMGetCoordinateDM(dm, &dmCoord);CHKERRQ(ierr);
613   ierr = SNESCreate(PETSC_COMM_SELF, &snes);CHKERRQ(ierr);
614   ierr = SNESSetOptionsPrefix(snes, "hex_interp_");CHKERRQ(ierr);
615   ierr = VecCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
616   ierr = VecSetSizes(r, 3, 3);CHKERRQ(ierr);
617   ierr = VecSetType(r,dm->vectype);CHKERRQ(ierr);
618   ierr = VecDuplicate(r, &ref);CHKERRQ(ierr);
619   ierr = VecDuplicate(r, &real);CHKERRQ(ierr);
620   ierr = MatCreate(PETSC_COMM_SELF, &J);CHKERRQ(ierr);
621   ierr = MatSetSizes(J, 3, 3, 3, 3);CHKERRQ(ierr);
622   ierr = MatSetType(J, MATSEQDENSE);CHKERRQ(ierr);
623   ierr = MatSetUp(J);CHKERRQ(ierr);
624   ierr = SNESSetFunction(snes, r, HexMap_Private, NULL);CHKERRQ(ierr);
625   ierr = SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);CHKERRQ(ierr);
626   ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
627   ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
628   ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
629   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
630 
631   ierr = VecGetArray(ctx->coords, &coords);CHKERRQ(ierr);
632   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
633   for (p = 0; p < ctx->n; ++p) {
634     PetscScalar *x = NULL, *vertices = NULL;
635     PetscScalar *xi;
636     PetscReal    xir[3];
637     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;
638 
639     /* Can make this do all points at once */
640     ierr = DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
641     if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);
642     ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
643     if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);
644     ierr   = SNESSetFunction(snes, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
645     ierr   = SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
646     ierr   = VecGetArray(real, &xi);CHKERRQ(ierr);
647     xi[0]  = coords[p*ctx->dim+0];
648     xi[1]  = coords[p*ctx->dim+1];
649     xi[2]  = coords[p*ctx->dim+2];
650     ierr   = VecRestoreArray(real, &xi);CHKERRQ(ierr);
651     ierr   = SNESSolve(snes, real, ref);CHKERRQ(ierr);
652     ierr   = VecGetArray(ref, &xi);CHKERRQ(ierr);
653     xir[0] = PetscRealPart(xi[0]);
654     xir[1] = PetscRealPart(xi[1]);
655     xir[2] = PetscRealPart(xi[2]);
656     for (comp = 0; comp < ctx->dof; ++comp) {
657       a[p*ctx->dof+comp] =
658         x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
659         x[3*ctx->dof+comp]*    xir[0]*(1-xir[1])*(1-xir[2]) +
660         x[2*ctx->dof+comp]*    xir[0]*    xir[1]*(1-xir[2]) +
661         x[1*ctx->dof+comp]*(1-xir[0])*    xir[1]*(1-xir[2]) +
662         x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*   xir[2] +
663         x[5*ctx->dof+comp]*    xir[0]*(1-xir[1])*   xir[2] +
664         x[6*ctx->dof+comp]*    xir[0]*    xir[1]*   xir[2] +
665         x[7*ctx->dof+comp]*(1-xir[0])*    xir[1]*   xir[2];
666     }
667     ierr = VecRestoreArray(ref, &xi);CHKERRQ(ierr);
668     ierr = DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
669     ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
670   }
671   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
672   ierr = VecRestoreArray(ctx->coords, &coords);CHKERRQ(ierr);
673 
674   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
675   ierr = VecDestroy(&r);CHKERRQ(ierr);
676   ierr = VecDestroy(&ref);CHKERRQ(ierr);
677   ierr = VecDestroy(&real);CHKERRQ(ierr);
678   ierr = MatDestroy(&J);CHKERRQ(ierr);
679   PetscFunctionReturn(0);
680 }
681 
682 #undef __FUNCT__
683 #define __FUNCT__ "DMInterpolationEvaluate"
684 /*
685   Input Parameters:
686 + ctx - The DMInterpolationInfo context
687 . dm  - The DM
688 - x   - The local vector containing the field to be interpolated
689 
690   Output Parameters:
691 . v   - The vector containing the interpolated values
692 */
693 PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
694 {
695   PetscInt       dim, coneSize, n;
696   PetscErrorCode ierr;
697 
698   PetscFunctionBegin;
699   PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
700   PetscValidHeaderSpecific(x, VEC_CLASSID, 3);
701   PetscValidHeaderSpecific(v, VEC_CLASSID, 4);
702   ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
703   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);
704   if (n) {
705     ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
706     ierr = DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);CHKERRQ(ierr);
707     if (dim == 2) {
708       if (coneSize == 3) {
709         ierr = DMInterpolate_Triangle_Private(ctx, dm, x, v);CHKERRQ(ierr);
710       } else if (coneSize == 4) {
711         ierr = DMInterpolate_Quad_Private(ctx, dm, x, v);CHKERRQ(ierr);
712       } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
713     } else if (dim == 3) {
714       if (coneSize == 4) {
715         ierr = DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);CHKERRQ(ierr);
716       } else {
717         ierr = DMInterpolate_Hex_Private(ctx, dm, x, v);CHKERRQ(ierr);
718       }
719     } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
720   }
721   PetscFunctionReturn(0);
722 }
723 
724 #undef __FUNCT__
725 #define __FUNCT__ "DMInterpolationDestroy"
726 PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
727 {
728   PetscErrorCode ierr;
729 
730   PetscFunctionBegin;
731   PetscValidPointer(ctx, 2);
732   ierr = VecDestroy(&(*ctx)->coords);CHKERRQ(ierr);
733   ierr = PetscFree((*ctx)->points);CHKERRQ(ierr);
734   ierr = PetscFree((*ctx)->cells);CHKERRQ(ierr);
735   ierr = PetscFree(*ctx);CHKERRQ(ierr);
736   *ctx = NULL;
737   PetscFunctionReturn(0);
738 }
739 
740 #undef __FUNCT__
741 #define __FUNCT__ "SNESMonitorFields"
742 /*@C
743   SNESMonitorFields - Monitors the residual for each field separately
744 
745   Collective on SNES
746 
747   Input Parameters:
748 + snes   - the SNES context
749 . its    - iteration number
750 . fgnorm - 2-norm of residual
751 - dummy  - unused context
752 
753   Notes:
754   This routine prints the residual norm at each iteration.
755 
756   Level: intermediate
757 
758 .keywords: SNES, nonlinear, default, monitor, norm
759 .seealso: SNESMonitorSet(), SNESMonitorDefault()
760 @*/
761 PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
762 {
763   PetscViewer        viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) snes));
764   Vec                res;
765   DM                 dm;
766   PetscSection       s;
767   const PetscScalar *r;
768   PetscReal         *lnorms, *norms;
769   PetscInt           numFields, f, pStart, pEnd, p;
770   PetscErrorCode     ierr;
771 
772   PetscFunctionBegin;
773   ierr = SNESGetFunction(snes, &res, 0, 0);CHKERRQ(ierr);
774   ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);
775   ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr);
776   ierr = PetscSectionGetNumFields(s, &numFields);CHKERRQ(ierr);
777   ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
778   ierr = PetscCalloc2(numFields, &lnorms, numFields, &norms);CHKERRQ(ierr);
779   ierr = VecGetArrayRead(res, &r);CHKERRQ(ierr);
780   for (p = pStart; p < pEnd; ++p) {
781     for (f = 0; f < numFields; ++f) {
782       PetscInt fdof, foff, d;
783 
784       ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr);
785       ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr);
786       for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
787     }
788   }
789   ierr = VecRestoreArrayRead(res, &r);CHKERRQ(ierr);
790   ierr = MPI_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr);
791   ierr = PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);CHKERRQ(ierr);
792   ierr = PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);CHKERRQ(ierr);
793   for (f = 0; f < numFields; ++f) {
794     if (f > 0) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);}
795     ierr = PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));CHKERRQ(ierr);
796   }
797   ierr = PetscViewerASCIIPrintf(viewer, "]\n");CHKERRQ(ierr);
798   ierr = PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);CHKERRQ(ierr);
799   ierr = PetscFree2(lnorms, norms);CHKERRQ(ierr);
800   PetscFunctionReturn(0);
801 }
802 
803 /********************* Residual Computation **************************/
804 
805 #undef __FUNCT__
806 #define __FUNCT__ "DMPlexSNESGetGeometryFEM"
807 /*@
808   DMPlexSNESGetGeometryFEM - Return precomputed geometric data
809 
810   Input Parameter:
811 . dm - The DM
812 
813   Output Parameters:
814 . cellgeom - The values precomputed from cell geometry
815 
816   Level: developer
817 
818 .seealso: DMPlexSNESSetFunctionLocal()
819 @*/
820 PetscErrorCode DMPlexSNESGetGeometryFEM(DM dm, Vec *cellgeom)
821 {
822   DMSNES         dmsnes;
823   PetscObject    obj;
824   PetscErrorCode ierr;
825 
826   PetscFunctionBegin;
827   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
828   ierr = DMGetDMSNES(dm, &dmsnes);CHKERRQ(ierr);
829   ierr = PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", &obj);CHKERRQ(ierr);
830   if (!obj) {
831     Vec cellgeom;
832 
833     ierr = DMPlexComputeGeometryFEM(dm, &cellgeom);CHKERRQ(ierr);
834     ierr = PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", (PetscObject) cellgeom);CHKERRQ(ierr);
835     ierr = VecDestroy(&cellgeom);CHKERRQ(ierr);
836   }
837   if (cellgeom) {PetscValidPointer(cellgeom, 3); ierr = PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", (PetscObject *) cellgeom);CHKERRQ(ierr);}
838   PetscFunctionReturn(0);
839 }
840 
841 #undef __FUNCT__
842 #define __FUNCT__ "DMPlexComputeResidual_Internal"
843 PetscErrorCode DMPlexComputeResidual_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
844 {
845   PetscErrorCode    ierr;
846 
847   PetscFunctionBegin;
848   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
849   /* FEM+FVM */
850   /* Get sizes from dm and dmAux */
851   /* ONCE: Compute geometric data (includes gradient), stash in mesh */
852   /*   Use named vectors */
853   /* Limit cell gradients */
854   /* Handle boundary values */
855   /* Loop over domain */
856   /*   Extract geometry and coefficients */
857   /* Loop over fields */
858   /*   Riemann solve over faces             (need fields at face centroids) */
859   /*   Integrate FE residual to get elemVec (need fields at quadrature points) */
860   /* Loop over domain */
861   /*   Add elemVec to locX */
862   /*   Accumulate fluxes to cells */
863 
864   /* FEM */
865   /* Get sizes from dm and dmAux */
866   /* Handle boundary values */
867   /* Loop over domain */
868   /*   Calculate geometry */
869   /*   Extract coefficients */
870   /* Loop over fields */
871   /*   Set tiling for FE*/
872   /*   Integrate FE residual to get elemVec */
873   /*     Loop over subdomain */
874   /*       Loop over quad points */
875   /*         Transform coords to real space */
876   /*         Evaluate field and aux fields at point */
877   /*         Evaluate residual at point */
878   /*         Transform residual to real space */
879   /*       Add residual to elemVec */
880   /* Loop over domain */
881   /*   Add elemVec to locX */
882 
883   /* FVM */
884   /* Compute geometric data */
885   /* If using gradients */
886   /*   Compute gradient data */
887   /*   Loop over domain faces */
888   /*     Count computational faces */
889   /*     Reconstruct cell gradient */
890   /*   Loop over domain cells */
891   /*     Limit cell gradients */
892   /* Handle boundary values */
893   /* Loop over domain faces */
894   /*   Read out field, centroid, normal, volume for each side of face */
895   /* Riemann solve over faces */
896   /* Loop over domain faces */
897   /*   Accumulate fluxes to cells */
898   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
899   PetscFunctionReturn(0);
900 }
901 
902 #undef __FUNCT__
903 #define __FUNCT__ "DMPlexComputeResidualFEM_Internal"
904 PetscErrorCode DMPlexComputeResidualFEM_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
905 {
906   DM_Plex          *mesh  = (DM_Plex *) dm->data;
907   const char       *name  = "Residual";
908   DM                dmAux;
909   DMLabel           depth;
910   Vec               A, cellgeom;
911   PetscDS           prob, probAux = NULL;
912   PetscQuadrature   q;
913   PetscSection      section, sectionAux;
914   PetscFECellGeom  *cgeom;
915   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
916   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c, numBd, bd;
917   PetscInt          totDim, totDimBd, totDimAux;
918   PetscErrorCode    ierr;
919 
920   PetscFunctionBegin;
921   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
922   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
923   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
924   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
925   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
926   ierr = PetscDSGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr);
927   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
928   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
929   numCells = cEnd - cStart;
930   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
931   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
932   if (dmAux) {
933     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
934     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
935     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
936   }
937   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
938   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
939   ierr = PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);CHKERRQ(ierr);
940   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
941   ierr = DMPlexSNESGetGeometryFEM(dm, &cellgeom);CHKERRQ(ierr);
942   ierr = VecGetArray(cellgeom, (PetscScalar **) &cgeom);CHKERRQ(ierr);
943   for (c = cStart; c < cEnd; ++c) {
944     PetscScalar *x = NULL, *x_t = NULL;
945     PetscInt     i;
946 
947     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
948     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
949     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
950     if (X_t) {
951       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
952       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
953       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
954     }
955     if (dmAux) {
956       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
957       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
958       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
959     }
960   }
961   for (f = 0; f < Nf; ++f) {
962     PetscFE  fe;
963     PetscInt numQuadPoints, Nb;
964     /* Conforming batches */
965     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
966     /* Remainder */
967     PetscInt Nr, offset;
968 
969     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
970     ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
971     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
972     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
973     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
974     blockSize = Nb*numQuadPoints;
975     batchSize = numBlocks * blockSize;
976     ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
977     numChunks = numCells / (numBatches*batchSize);
978     Ne        = numChunks*numBatches*batchSize;
979     Nr        = numCells % (numBatches*batchSize);
980     offset    = numCells - Nr;
981     ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVec);CHKERRQ(ierr);
982     ierr = PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);CHKERRQ(ierr);
983   }
984   for (c = cStart; c < cEnd; ++c) {
985     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, totDim, &elemVec[c*totDim]);CHKERRQ(ierr);}
986     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);CHKERRQ(ierr);
987   }
988   ierr = VecRestoreArray(cellgeom, (PetscScalar **) &cgeom);CHKERRQ(ierr);
989   ierr = PetscFree3(u,u_t,elemVec);CHKERRQ(ierr);
990   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
991   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
992   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
993   for (bd = 0; bd < numBd; ++bd) {
994     const char     *bdLabel;
995     DMLabel         label;
996     IS              pointIS;
997     const PetscInt *points;
998     const PetscInt *values;
999     PetscInt        field, numValues, numPoints, p, dep, numFaces;
1000     PetscBool       isEssential;
1001 
1002     ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr);
1003     if (isEssential) continue;
1004     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1005     ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);
1006     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
1007     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
1008     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
1009     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1010       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1011       if (dep == dim-1) ++numFaces;
1012     }
1013     ierr = PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd,&elemVec);CHKERRQ(ierr);
1014     if (X_t) {ierr = PetscMalloc1(numFaces*totDimBd,&u_t);CHKERRQ(ierr);}
1015     for (p = 0, f = 0; p < numPoints; ++p) {
1016       const PetscInt point = points[p];
1017       PetscScalar   *x     = NULL;
1018       PetscInt       i;
1019 
1020       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1021       if (dep != dim-1) continue;
1022       ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);CHKERRQ(ierr);
1023       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
1024       if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
1025       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1026       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1027       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1028       if (X_t) {
1029         ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1030         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1031         ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1032       }
1033       ++f;
1034     }
1035     for (f = 0; f < Nf; ++f) {
1036       PetscFE  fe;
1037       PetscInt numQuadPoints, Nb;
1038       /* Conforming batches */
1039       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1040       /* Remainder */
1041       PetscInt Nr, offset;
1042 
1043       ierr = PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1044       ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1045       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1046       ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1047       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1048       blockSize = Nb*numQuadPoints;
1049       batchSize = numBlocks * blockSize;
1050       ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1051       numChunks = numFaces / (numBatches*batchSize);
1052       Ne        = numChunks*numBatches*batchSize;
1053       Nr        = numFaces % (numBatches*batchSize);
1054       offset    = numFaces - Nr;
1055       ierr = PetscFEIntegrateBdResidual(fe, prob, f, Ne, cgeom, u, u_t, NULL, NULL, elemVec);CHKERRQ(ierr);
1056       ierr = PetscFEIntegrateBdResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);CHKERRQ(ierr);
1057     }
1058     for (p = 0, f = 0; p < numPoints; ++p) {
1059       const PetscInt point = points[p];
1060 
1061       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1062       if (dep != dim-1) continue;
1063       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);CHKERRQ(ierr);}
1064       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*totDimBd], ADD_VALUES);CHKERRQ(ierr);
1065       ++f;
1066     }
1067     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1068     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1069     ierr = PetscFree3(u,cgeom,elemVec);CHKERRQ(ierr);
1070     if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);}
1071   }
1072   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
1073   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
1074   PetscFunctionReturn(0);
1075 }
1076 
1077 #undef __FUNCT__
1078 #define __FUNCT__ "DMPlexComputeResidualFEM_Check_Internal"
1079 static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
1080 {
1081   DM                dmCh, dmAux;
1082   Vec               A, cellgeom;
1083   PetscDS           prob, probCh, probAux = NULL;
1084   PetscQuadrature   q;
1085   PetscSection      section, sectionAux;
1086   PetscFECellGeom  *cgeom;
1087   PetscScalar      *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1088   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
1089   PetscInt          totDim, totDimAux, diffCell = 0;
1090   PetscErrorCode    ierr;
1091 
1092   PetscFunctionBegin;
1093   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1094   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1095   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1096   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1097   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1098   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1099   numCells = cEnd - cStart;
1100   ierr = PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);CHKERRQ(ierr);
1101   ierr = DMGetDS(dmCh, &probCh);CHKERRQ(ierr);
1102   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1103   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1104   if (dmAux) {
1105     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1106     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1107     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1108   }
1109   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
1110   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
1111   ierr = PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);CHKERRQ(ierr);
1112   ierr = PetscMalloc1(numCells*totDim,&elemVecCh);CHKERRQ(ierr);
1113   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1114   ierr = DMPlexSNESGetGeometryFEM(dm, &cellgeom);CHKERRQ(ierr);
1115   ierr = VecGetArray(cellgeom, (PetscScalar **) &cgeom);CHKERRQ(ierr);
1116   for (c = cStart; c < cEnd; ++c) {
1117     PetscScalar *x = NULL, *x_t = NULL;
1118     PetscInt     i;
1119 
1120     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1121     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1122     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1123     if (X_t) {
1124       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1125       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1126       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1127     }
1128     if (dmAux) {
1129       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1130       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1131       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1132     }
1133   }
1134   for (f = 0; f < Nf; ++f) {
1135     PetscFE  fe, feCh;
1136     PetscInt numQuadPoints, Nb;
1137     /* Conforming batches */
1138     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1139     /* Remainder */
1140     PetscInt Nr, offset;
1141 
1142     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1143     ierr = PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);CHKERRQ(ierr);
1144     ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1145     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1146     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1147     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1148     blockSize = Nb*numQuadPoints;
1149     batchSize = numBlocks * blockSize;
1150     ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1151     numChunks = numCells / (numBatches*batchSize);
1152     Ne        = numChunks*numBatches*batchSize;
1153     Nr        = numCells % (numBatches*batchSize);
1154     offset    = numCells - Nr;
1155     ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVec);CHKERRQ(ierr);
1156     ierr = PetscFEIntegrateResidual(feCh, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVecCh);CHKERRQ(ierr);
1157     ierr = PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);CHKERRQ(ierr);
1158     ierr = PetscFEIntegrateResidual(feCh, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVecCh[offset*totDim]);CHKERRQ(ierr);
1159   }
1160   for (c = cStart; c < cEnd; ++c) {
1161     PetscBool diff = PETSC_FALSE;
1162     PetscInt  d;
1163 
1164     for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1165     if (diff) {
1166       ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);CHKERRQ(ierr);
1167       ierr = DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);CHKERRQ(ierr);
1168       ierr = DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);CHKERRQ(ierr);
1169       ++diffCell;
1170     }
1171     if (diffCell > 9) break;
1172     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);CHKERRQ(ierr);
1173   }
1174   ierr = VecRestoreArray(cellgeom, (PetscScalar **) &cgeom);CHKERRQ(ierr);
1175   ierr = PetscFree3(u,u_t,elemVec);CHKERRQ(ierr);
1176   ierr = PetscFree(elemVecCh);CHKERRQ(ierr);
1177   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1178   PetscFunctionReturn(0);
1179 }
1180 
1181 #undef __FUNCT__
1182 #define __FUNCT__ "DMPlexSNESComputeResidualFEM"
1183 /*@
1184   DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
1185 
1186   Input Parameters:
1187 + dm - The mesh
1188 . X  - Local solution
1189 - user - The user context
1190 
1191   Output Parameter:
1192 . F  - Local output vector
1193 
1194   Level: developer
1195 
1196 .seealso: DMPlexComputeJacobianActionFEM()
1197 @*/
1198 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1199 {
1200   PetscObject    check;
1201   PetscErrorCode ierr;
1202 
1203   PetscFunctionBegin;
1204   /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
1205   ierr = PetscObjectQuery((PetscObject) dm, "dmCh", &check);CHKERRQ(ierr);
1206   if (check) {ierr = DMPlexComputeResidualFEM_Check_Internal(dm, X, NULL, F, user);CHKERRQ(ierr);}
1207   else       {ierr = DMPlexComputeResidualFEM_Internal(dm, X, NULL, F, user);CHKERRQ(ierr);}
1208   PetscFunctionReturn(0);
1209 }
1210 
1211 #undef __FUNCT__
1212 #define __FUNCT__ "DMPlexComputeJacobianFEM_Internal"
1213 PetscErrorCode DMPlexComputeJacobianFEM_Internal(DM dm, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user)
1214 {
1215   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1216   const char       *name  = "Jacobian";
1217   DM                dmAux;
1218   DMLabel           depth;
1219   Vec               A, cellgeom;
1220   PetscDS           prob, probAux = NULL;
1221   PetscQuadrature   quad;
1222   PetscSection      section, globalSection, sectionAux;
1223   PetscFECellGeom  *cgeom;
1224   PetscScalar      *elemMat, *u, *u_t, *a = NULL;
1225   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1226   PetscInt          totDim, totDimBd, totDimAux, numBd, bd;
1227   PetscBool         isShell;
1228   PetscErrorCode    ierr;
1229 
1230   PetscFunctionBegin;
1231   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1232   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1233   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1234   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1235   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1236   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1237   ierr = PetscDSGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr);
1238   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1239   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1240   numCells = cEnd - cStart;
1241   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1242   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1243   if (dmAux) {
1244     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1245     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1246     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1247   }
1248   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
1249   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1250   ierr = PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat);CHKERRQ(ierr);
1251   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1252   ierr = DMPlexSNESGetGeometryFEM(dm, &cellgeom);CHKERRQ(ierr);
1253   ierr = VecGetArray(cellgeom, (PetscScalar **) &cgeom);CHKERRQ(ierr);
1254   for (c = cStart; c < cEnd; ++c) {
1255     PetscScalar *x = NULL,  *x_t = NULL;
1256     PetscInt     i;
1257 
1258     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1259     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1260     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1261     if (X_t) {
1262       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1263       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1264       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1265     }
1266     if (dmAux) {
1267       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1268       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1269       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1270     }
1271   }
1272   ierr = PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));CHKERRQ(ierr);
1273   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1274     PetscFE  fe;
1275     PetscInt numQuadPoints, Nb;
1276     /* Conforming batches */
1277     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1278     /* Remainder */
1279     PetscInt Nr, offset;
1280 
1281     ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1282     ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1283     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1284     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1285     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1286     blockSize = Nb*numQuadPoints;
1287     batchSize = numBlocks * blockSize;
1288     ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1289     numChunks = numCells / (numBatches*batchSize);
1290     Ne        = numChunks*numBatches*batchSize;
1291     Nr        = numCells % (numBatches*batchSize);
1292     offset    = numCells - Nr;
1293     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1294       ierr = PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMat);CHKERRQ(ierr);
1295       ierr = 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]);CHKERRQ(ierr);
1296     }
1297   }
1298   for (c = cStart; c < cEnd; ++c) {
1299     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[c*totDim*totDim]);CHKERRQ(ierr);}
1300     ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*totDim*totDim], ADD_VALUES);CHKERRQ(ierr);
1301   }
1302   ierr = VecGetArray(cellgeom, (PetscScalar **) &cgeom);CHKERRQ(ierr);
1303   ierr = PetscFree3(u,u_t,elemMat);CHKERRQ(ierr);
1304   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1305   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
1306   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
1307   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
1308   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
1309   for (bd = 0; bd < numBd; ++bd) {
1310     const char     *bdLabel;
1311     DMLabel         label;
1312     IS              pointIS;
1313     const PetscInt *points;
1314     const PetscInt *values;
1315     PetscInt        field, numValues, numPoints, p, dep, numFaces;
1316     PetscBool       isEssential;
1317 
1318     ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr);
1319     if (isEssential) continue;
1320     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1321     ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);
1322     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
1323     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
1324     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
1325     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1326       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1327       if (dep == dim-1) ++numFaces;
1328     }
1329     ierr = PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd*totDimBd,&elemMat);CHKERRQ(ierr);
1330     if (X_t) {ierr = PetscMalloc1(numFaces*totDimBd,&u_t);CHKERRQ(ierr);}
1331     for (p = 0, f = 0; p < numPoints; ++p) {
1332       const PetscInt point = points[p];
1333       PetscScalar   *x     = NULL;
1334       PetscInt       i;
1335 
1336       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1337       if (dep != dim-1) continue;
1338       ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);CHKERRQ(ierr);
1339       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
1340       if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
1341       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1342       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1343       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1344       if (X_t) {
1345         ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1346         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1347         ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1348       }
1349       ++f;
1350     }
1351     ierr = PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));CHKERRQ(ierr);
1352     for (fieldI = 0; fieldI < Nf; ++fieldI) {
1353       PetscFE  fe;
1354       PetscInt numQuadPoints, Nb;
1355       /* Conforming batches */
1356       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1357       /* Remainder */
1358       PetscInt Nr, offset;
1359 
1360       ierr = PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1361       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1362       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1363       ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1364       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1365       blockSize = Nb*numQuadPoints;
1366       batchSize = numBlocks * blockSize;
1367       ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1368       numChunks = numFaces / (numBatches*batchSize);
1369       Ne        = numChunks*numBatches*batchSize;
1370       Nr        = numFaces % (numBatches*batchSize);
1371       offset    = numFaces - Nr;
1372       for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1373         ierr = PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, NULL, NULL, elemMat);CHKERRQ(ierr);
1374         ierr = PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);CHKERRQ(ierr);
1375       }
1376     }
1377     for (p = 0, f = 0; p < numPoints; ++p) {
1378       const PetscInt point = points[p];
1379 
1380       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1381       if (dep != dim-1) continue;
1382       if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);CHKERRQ(ierr);}
1383       ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);CHKERRQ(ierr);
1384       ++f;
1385     }
1386     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1387     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1388     ierr = PetscFree3(u,cgeom,elemMat);CHKERRQ(ierr);
1389     if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);}
1390   }
1391   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1392   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1393   if (mesh->printFEM) {
1394     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1395     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
1396     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1397   }
1398   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1399   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
1400   if (isShell) {
1401     JacActionCtx *jctx;
1402 
1403     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1404     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
1405   }
1406   PetscFunctionReturn(0);
1407 }
1408 
1409 #undef __FUNCT__
1410 #define __FUNCT__ "DMPlexSNESComputeJacobianFEM"
1411 /*@
1412   DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1413 
1414   Input Parameters:
1415 + dm - The mesh
1416 . X  - Local input vector
1417 - user - The user context
1418 
1419   Output Parameter:
1420 . Jac  - Jacobian matrix
1421 
1422   Note:
1423   The first member of the user context must be an FEMContext.
1424 
1425   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1426   like a GPU, or vectorize on a multicore machine.
1427 
1428   Level: developer
1429 
1430 .seealso: FormFunctionLocal()
1431 @*/
1432 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1433 {
1434   PetscErrorCode ierr;
1435 
1436   PetscFunctionBegin;
1437   ierr = DMPlexComputeJacobianFEM_Internal(dm, X, NULL, Jac, JacP, user);CHKERRQ(ierr);
1438   PetscFunctionReturn(0);
1439 }
1440