xref: /petsc/src/snes/utils/dmplexsnes.c (revision 24cdb843c62e1437a3052e6f6c6af2ebc2f9cd7c) 
1 #include <petsc-private/dmpleximpl.h>   /*I "petscdmplex.h" I*/
2 #include <petscsnes.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__ "DMPlexComputeResidual_Internal"
807 PetscErrorCode DMPlexComputeResidual_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
808 {
809   PetscErrorCode    ierr;
810 
811   PetscFunctionBegin;
812   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
813   /* FEM+FVM */
814   /* Get sizes from dm and dmAux */
815   /* ONCE: Compute geometric data (includes gradient), stash in mesh */
816   /*   Use named vectors */
817   /* Limit cell gradients */
818   /* Handle boundary values */
819   /* Loop over domain */
820   /*   Extract geometry and coefficients */
821   /* Loop over fields */
822   /*   Riemann solve over faces             (need fields at face centroids) */
823   /*   Integrate FE residual to get elemVec (need fields at quadrature points) */
824   /* Loop over domain */
825   /*   Add elemVec to locX */
826   /*   Accumulate fluxes to cells */
827 
828   /* FEM */
829   /* Get sizes from dm and dmAux */
830   /* Handle boundary values */
831   /* Loop over domain */
832   /*   Calculate geometry */
833   /*   Extract coefficients */
834   /* Loop over fields */
835   /*   Set tiling for FE*/
836   /*   Integrate FE residual to get elemVec */
837   /*     Loop over subdomain */
838   /*       Loop over quad points */
839   /*         Transform coords to real space */
840   /*         Evaluate field and aux fields at point */
841   /*         Evaluate residual at point */
842   /*         Transform residual to real space */
843   /*       Add residual to elemVec */
844   /* Loop over domain */
845   /*   Add elemVec to locX */
846 
847   /* FVM */
848   /* Compute geometric data */
849   /* If using gradients */
850   /*   Compute gradient data */
851   /*   Loop over domain faces */
852   /*     Count computational faces */
853   /*     Reconstruct cell gradient */
854   /*   Loop over domain cells */
855   /*     Limit cell gradients */
856   /* Handle boundary values */
857   /* Loop over domain faces */
858   /*   Read out field, centroid, normal, volume for each side of face */
859   /* Riemann solve over faces */
860   /* Loop over domain faces */
861   /*   Accumulate fluxes to cells */
862   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
863   PetscFunctionReturn(0);
864 }
865 
866 #undef __FUNCT__
867 #define __FUNCT__ "DMPlexComputeResidualFEM_Internal"
868 PetscErrorCode DMPlexComputeResidualFEM_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
869 {
870   DM_Plex          *mesh  = (DM_Plex *) dm->data;
871   const char       *name  = "Residual";
872   DM                dmAux;
873   DMLabel           depth;
874   Vec               A;
875   PetscDS           prob, probAux = NULL;
876   PetscQuadrature   q;
877   PetscCellGeometry geom;
878   PetscSection      section, sectionAux;
879   PetscReal        *v0, *J, *invJ, *detJ;
880   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
881   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c, numBd, bd;
882   PetscInt          totDim, totDimBd, totDimAux;
883   PetscErrorCode    ierr;
884 
885   PetscFunctionBegin;
886   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
887   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
888   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
889   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
890   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
891   ierr = PetscDSGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr);
892   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
893   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
894   numCells = cEnd - cStart;
895   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
896   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
897   if (dmAux) {
898     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
899     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
900     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
901   }
902   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
903   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
904   ierr = PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim,&elemVec);CHKERRQ(ierr);
905   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
906   for (c = cStart; c < cEnd; ++c) {
907     PetscScalar *x = NULL, *x_t = NULL;
908     PetscInt     i;
909 
910     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
911     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
912     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
913     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
914     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
915     if (X_t) {
916       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
917       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
918       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
919     }
920     if (dmAux) {
921       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
922       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
923       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
924     }
925   }
926   for (f = 0; f < Nf; ++f) {
927     PetscFE  fe;
928     PetscInt numQuadPoints, Nb;
929     /* Conforming batches */
930     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
931     /* Remainder */
932     PetscInt Nr, offset;
933 
934     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
935     ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
936     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
937     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
938     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
939     blockSize = Nb*numQuadPoints;
940     batchSize = numBlocks * blockSize;
941     ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
942     numChunks = numCells / (numBatches*batchSize);
943     Ne        = numChunks*numBatches*batchSize;
944     Nr        = numCells % (numBatches*batchSize);
945     offset    = numCells - Nr;
946     geom.v0   = v0;
947     geom.J    = J;
948     geom.invJ = invJ;
949     geom.detJ = detJ;
950     ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, geom, u, u_t, probAux, a, elemVec);CHKERRQ(ierr);
951     geom.v0   = &v0[offset*dim];
952     geom.J    = &J[offset*dim*dim];
953     geom.invJ = &invJ[offset*dim*dim];
954     geom.detJ = &detJ[offset];
955     ierr = PetscFEIntegrateResidual(fe, prob, f, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);CHKERRQ(ierr);
956   }
957   for (c = cStart; c < cEnd; ++c) {
958     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, totDim, &elemVec[c*totDim]);CHKERRQ(ierr);}
959     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);CHKERRQ(ierr);
960   }
961   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
962   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
963   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
964   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
965   for (bd = 0; bd < numBd; ++bd) {
966     const char     *bdLabel;
967     DMLabel         label;
968     IS              pointIS;
969     const PetscInt *points;
970     const PetscInt *values;
971     PetscReal      *n;
972     PetscInt        field, numValues, numPoints, p, dep, numFaces;
973     PetscBool       isEssential;
974 
975     ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr);
976     if (isEssential) continue;
977     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
978     ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);
979     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
980     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
981     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
982     for (p = 0, numFaces = 0; p < numPoints; ++p) {
983       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
984       if (dep == dim-1) ++numFaces;
985     }
986     ierr = PetscMalloc7(numFaces*totDimBd,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*totDimBd,&elemVec);CHKERRQ(ierr);
987     if (X_t) {ierr = PetscMalloc1(numFaces*totDimBd,&u_t);CHKERRQ(ierr);}
988     for (p = 0, f = 0; p < numPoints; ++p) {
989       const PetscInt point = points[p];
990       PetscScalar   *x     = NULL;
991       PetscInt       i;
992 
993       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
994       if (dep != dim-1) continue;
995       ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
996       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
997       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
998       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
999       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1000       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1001       if (X_t) {
1002         ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1003         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1004         ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1005       }
1006       ++f;
1007     }
1008     for (f = 0; f < Nf; ++f) {
1009       PetscFE  fe;
1010       PetscInt numQuadPoints, Nb;
1011       /* Conforming batches */
1012       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1013       /* Remainder */
1014       PetscInt Nr, offset;
1015 
1016       ierr = PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1017       ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1018       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1019       ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1020       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1021       blockSize = Nb*numQuadPoints;
1022       batchSize = numBlocks * blockSize;
1023       ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1024       numChunks = numFaces / (numBatches*batchSize);
1025       Ne        = numChunks*numBatches*batchSize;
1026       Nr        = numFaces % (numBatches*batchSize);
1027       offset    = numFaces - Nr;
1028       geom.v0   = v0;
1029       geom.n    = n;
1030       geom.J    = J;
1031       geom.invJ = invJ;
1032       geom.detJ = detJ;
1033       ierr = PetscFEIntegrateBdResidual(fe, prob, f, Ne, geom, u, u_t, NULL, NULL, elemVec);CHKERRQ(ierr);
1034       geom.v0   = &v0[offset*dim];
1035       geom.n    = &n[offset*dim];
1036       geom.J    = &J[offset*dim*dim];
1037       geom.invJ = &invJ[offset*dim*dim];
1038       geom.detJ = &detJ[offset];
1039       ierr = PetscFEIntegrateBdResidual(fe, prob, f, Nr, geom, &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);CHKERRQ(ierr);
1040     }
1041     for (p = 0, f = 0; p < numPoints; ++p) {
1042       const PetscInt point = points[p];
1043 
1044       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1045       if (dep != dim-1) continue;
1046       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);CHKERRQ(ierr);}
1047       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*totDimBd], ADD_VALUES);CHKERRQ(ierr);
1048       ++f;
1049     }
1050     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1051     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1052     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1053     if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);}
1054   }
1055   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
1056   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
1057   PetscFunctionReturn(0);
1058 }
1059 
1060 #undef __FUNCT__
1061 #define __FUNCT__ "DMPlexComputeResidualFEM_Check_Internal"
1062 static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
1063 {
1064   DM                dmCh, dmAux;
1065   Vec               A;
1066   PetscDS           prob, probCh, probAux = NULL;
1067   PetscQuadrature   q;
1068   PetscCellGeometry geom;
1069   PetscSection      section, sectionAux;
1070   PetscReal        *v0, *J, *invJ, *detJ;
1071   PetscScalar      *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1072   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
1073   PetscInt          totDim, totDimAux, diffCell = 0;
1074   PetscErrorCode    ierr;
1075 
1076   PetscFunctionBegin;
1077   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1078   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1079   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1080   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1081   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1082   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1083   numCells = cEnd - cStart;
1084   ierr = PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);CHKERRQ(ierr);
1085   ierr = DMGetDS(dmCh, &probCh);CHKERRQ(ierr);
1086   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1087   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1088   if (dmAux) {
1089     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1090     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1091     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1092   }
1093   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
1094   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
1095   ierr = PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim,&elemVec);CHKERRQ(ierr);
1096   ierr = PetscMalloc1(numCells*totDim,&elemVecCh);CHKERRQ(ierr);
1097   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1098   for (c = cStart; c < cEnd; ++c) {
1099     PetscScalar *x = NULL, *x_t = NULL;
1100     PetscInt     i;
1101 
1102     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1103     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1104     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1105     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1106     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1107     if (X_t) {
1108       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1109       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1110       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1111     }
1112     if (dmAux) {
1113       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1114       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1115       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1116     }
1117   }
1118   for (f = 0; f < Nf; ++f) {
1119     PetscFE  fe, feCh;
1120     PetscInt numQuadPoints, Nb;
1121     /* Conforming batches */
1122     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1123     /* Remainder */
1124     PetscInt Nr, offset;
1125 
1126     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1127     ierr = PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);CHKERRQ(ierr);
1128     ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1129     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1130     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1131     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1132     blockSize = Nb*numQuadPoints;
1133     batchSize = numBlocks * blockSize;
1134     ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1135     numChunks = numCells / (numBatches*batchSize);
1136     Ne        = numChunks*numBatches*batchSize;
1137     Nr        = numCells % (numBatches*batchSize);
1138     offset    = numCells - Nr;
1139     geom.v0   = v0;
1140     geom.J    = J;
1141     geom.invJ = invJ;
1142     geom.detJ = detJ;
1143     ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, geom, u, u_t, probAux, a, elemVec);CHKERRQ(ierr);
1144     ierr = PetscFEIntegrateResidual(feCh, prob, f, Ne, geom, u, u_t, probAux, a, elemVecCh);CHKERRQ(ierr);
1145     geom.v0   = &v0[offset*dim];
1146     geom.J    = &J[offset*dim*dim];
1147     geom.invJ = &invJ[offset*dim*dim];
1148     geom.detJ = &detJ[offset];
1149     ierr = PetscFEIntegrateResidual(fe, prob, f, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);CHKERRQ(ierr);
1150     ierr = PetscFEIntegrateResidual(feCh, prob, f, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVecCh[offset*totDim]);CHKERRQ(ierr);
1151   }
1152   for (c = cStart; c < cEnd; ++c) {
1153     PetscBool diff = PETSC_FALSE;
1154     PetscInt  d;
1155 
1156     for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1157     if (diff) {
1158       ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);CHKERRQ(ierr);
1159       ierr = DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);CHKERRQ(ierr);
1160       ierr = DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);CHKERRQ(ierr);
1161       ++diffCell;
1162     }
1163     if (diffCell > 9) break;
1164     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);CHKERRQ(ierr);
1165   }
1166   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1167   ierr = PetscFree(elemVecCh);CHKERRQ(ierr);
1168   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1169   PetscFunctionReturn(0);
1170 }
1171 
1172 #undef __FUNCT__
1173 #define __FUNCT__ "DMPlexSNESComputeResidualFEM"
1174 /*@
1175   DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
1176 
1177   Input Parameters:
1178 + dm - The mesh
1179 . X  - Local solution
1180 - user - The user context
1181 
1182   Output Parameter:
1183 . F  - Local output vector
1184 
1185   Level: developer
1186 
1187 .seealso: DMPlexComputeJacobianActionFEM()
1188 @*/
1189 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1190 {
1191   PetscObject    check;
1192   PetscErrorCode ierr;
1193 
1194   PetscFunctionBegin;
1195   /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
1196   ierr = PetscObjectQuery((PetscObject) dm, "dmCh", &check);CHKERRQ(ierr);
1197   if (check) {ierr = DMPlexComputeResidualFEM_Check_Internal(dm, X, NULL, F, user);CHKERRQ(ierr);}
1198   else       {ierr = DMPlexComputeResidualFEM_Internal(dm, X, NULL, F, user);CHKERRQ(ierr);}
1199   PetscFunctionReturn(0);
1200 }
1201 
1202 #undef __FUNCT__
1203 #define __FUNCT__ "DMPlexComputeJacobianFEM_Internal"
1204 PetscErrorCode DMPlexComputeJacobianFEM_Internal(DM dm, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user)
1205 {
1206   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1207   const char       *name  = "Jacobian";
1208   DM                dmAux;
1209   DMLabel           depth;
1210   Vec               A;
1211   PetscDS           prob, probAux = NULL;
1212   PetscQuadrature   quad;
1213   PetscCellGeometry geom;
1214   PetscSection      section, globalSection, sectionAux;
1215   PetscReal        *v0, *J, *invJ, *detJ;
1216   PetscScalar      *elemMat, *u, *u_t, *a = NULL;
1217   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1218   PetscInt          totDim, totDimBd, totDimAux, numBd, bd;
1219   PetscBool         isShell;
1220   PetscErrorCode    ierr;
1221 
1222   PetscFunctionBegin;
1223   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1224   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1225   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1226   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1227   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1228   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1229   ierr = PetscDSGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr);
1230   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1231   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1232   numCells = cEnd - cStart;
1233   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1234   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1235   if (dmAux) {
1236     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1237     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1238     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1239   }
1240   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
1241   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1242   ierr = PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim*totDim,&elemMat);CHKERRQ(ierr);
1243   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1244   for (c = cStart; c < cEnd; ++c) {
1245     PetscScalar *x = NULL,  *x_t = NULL;
1246     PetscInt     i;
1247 
1248     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1249     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1250     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1251     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1252     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1253     if (X_t) {
1254       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1255       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1256       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1257     }
1258     if (dmAux) {
1259       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1260       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1261       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1262     }
1263   }
1264   ierr = PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));CHKERRQ(ierr);
1265   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1266     PetscFE  fe;
1267     PetscInt numQuadPoints, Nb;
1268     /* Conforming batches */
1269     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1270     /* Remainder */
1271     PetscInt Nr, offset;
1272 
1273     ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1274     ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1275     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1276     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1277     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1278     blockSize = Nb*numQuadPoints;
1279     batchSize = numBlocks * blockSize;
1280     ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1281     numChunks = numCells / (numBatches*batchSize);
1282     Ne        = numChunks*numBatches*batchSize;
1283     Nr        = numCells % (numBatches*batchSize);
1284     offset    = numCells - Nr;
1285     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1286       geom.v0   = v0;
1287       geom.J    = J;
1288       geom.invJ = invJ;
1289       geom.detJ = detJ;
1290       ierr = PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, probAux, a, elemMat);CHKERRQ(ierr);
1291       geom.v0   = &v0[offset*dim];
1292       geom.J    = &J[offset*dim*dim];
1293       geom.invJ = &invJ[offset*dim*dim];
1294       geom.detJ = &detJ[offset];
1295       ierr = PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Nr, geom, &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 = PetscFree7(u,u_t,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1303   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1304   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
1305   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
1306   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
1307   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
1308   for (bd = 0; bd < numBd; ++bd) {
1309     const char     *bdLabel;
1310     DMLabel         label;
1311     IS              pointIS;
1312     const PetscInt *points;
1313     const PetscInt *values;
1314     PetscReal      *n;
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 = PetscMalloc7(numFaces*totDimBd,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,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, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
1339       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
1340       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], 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         geom.v0   = v0;
1374         geom.n    = n;
1375         geom.J    = J;
1376         geom.invJ = invJ;
1377         geom.detJ = detJ;
1378         ierr = PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, NULL, NULL, elemMat);CHKERRQ(ierr);
1379         geom.v0   = &v0[offset*dim];
1380         geom.n    = &n[offset*dim];
1381         geom.J    = &J[offset*dim*dim];
1382         geom.invJ = &invJ[offset*dim*dim];
1383         geom.detJ = &detJ[offset];
1384         ierr = PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, geom, &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);CHKERRQ(ierr);
1385       }
1386     }
1387     for (p = 0, f = 0; p < numPoints; ++p) {
1388       const PetscInt point = points[p];
1389 
1390       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1391       if (dep != dim-1) continue;
1392       if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);CHKERRQ(ierr);}
1393       ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);CHKERRQ(ierr);
1394       ++f;
1395     }
1396     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1397     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1398     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1399     if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);}
1400   }
1401   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1402   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1403   if (mesh->printFEM) {
1404     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1405     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
1406     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1407   }
1408   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1409   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
1410   if (isShell) {
1411     JacActionCtx *jctx;
1412 
1413     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1414     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
1415   }
1416   PetscFunctionReturn(0);
1417 }
1418 
1419 #undef __FUNCT__
1420 #define __FUNCT__ "DMPlexSNESComputeJacobianFEM"
1421 /*@
1422   DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1423 
1424   Input Parameters:
1425 + dm - The mesh
1426 . X  - Local input vector
1427 - user - The user context
1428 
1429   Output Parameter:
1430 . Jac  - Jacobian matrix
1431 
1432   Note:
1433   The first member of the user context must be an FEMContext.
1434 
1435   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1436   like a GPU, or vectorize on a multicore machine.
1437 
1438   Level: developer
1439 
1440 .seealso: FormFunctionLocal()
1441 @*/
1442 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1443 {
1444   PetscErrorCode ierr;
1445 
1446   PetscFunctionBegin;
1447   ierr = DMPlexComputeJacobianFEM_Internal(dm, X, NULL, Jac, JacP, user);CHKERRQ(ierr);
1448   PetscFunctionReturn(0);
1449 }
1450