xref: /petsc/src/snes/utils/dmplexsnes.c (revision 7d4028c8da78af24e52f42b016cff0f57ae33c4d) 
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;
911   PetscDS           prob, probAux = NULL;
912   PetscQuadrature   q;
913   PetscCellGeometry geom;
914   PetscSection      section, sectionAux;
915   PetscReal        *v0, *J, *invJ, *detJ;
916   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
917   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c, numBd, bd;
918   PetscInt          totDim, totDimBd, totDimAux;
919   PetscErrorCode    ierr;
920 
921   PetscFunctionBegin;
922   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
923   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
924   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
925   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
926   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
927   ierr = PetscDSGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr);
928   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
929   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
930   numCells = cEnd - cStart;
931   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
932   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
933   if (dmAux) {
934     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
935     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
936     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
937   }
938   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
939   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
940   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);
941   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
942   for (c = cStart; c < cEnd; ++c) {
943     PetscScalar *x = NULL, *x_t = NULL;
944     PetscInt     i;
945 
946     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
947     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
948     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
949     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
950     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
951     if (X_t) {
952       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
953       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
954       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
955     }
956     if (dmAux) {
957       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
958       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
959       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
960     }
961   }
962   for (f = 0; f < Nf; ++f) {
963     PetscFE  fe;
964     PetscInt numQuadPoints, Nb;
965     /* Conforming batches */
966     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
967     /* Remainder */
968     PetscInt Nr, offset;
969 
970     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
971     ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
972     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
973     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
974     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
975     blockSize = Nb*numQuadPoints;
976     batchSize = numBlocks * blockSize;
977     ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
978     numChunks = numCells / (numBatches*batchSize);
979     Ne        = numChunks*numBatches*batchSize;
980     Nr        = numCells % (numBatches*batchSize);
981     offset    = numCells - Nr;
982     geom.v0   = v0;
983     geom.J    = J;
984     geom.invJ = invJ;
985     geom.detJ = detJ;
986     ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, geom, u, u_t, probAux, a, elemVec);CHKERRQ(ierr);
987     geom.v0   = &v0[offset*dim];
988     geom.J    = &J[offset*dim*dim];
989     geom.invJ = &invJ[offset*dim*dim];
990     geom.detJ = &detJ[offset];
991     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);
992   }
993   for (c = cStart; c < cEnd; ++c) {
994     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, totDim, &elemVec[c*totDim]);CHKERRQ(ierr);}
995     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);CHKERRQ(ierr);
996   }
997   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
998   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
999   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
1000   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
1001   for (bd = 0; bd < numBd; ++bd) {
1002     const char     *bdLabel;
1003     DMLabel         label;
1004     IS              pointIS;
1005     const PetscInt *points;
1006     const PetscInt *values;
1007     PetscReal      *n;
1008     PetscInt        field, numValues, numPoints, p, dep, numFaces;
1009     PetscBool       isEssential;
1010 
1011     ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr);
1012     if (isEssential) continue;
1013     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1014     ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);
1015     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
1016     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
1017     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
1018     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1019       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1020       if (dep == dim-1) ++numFaces;
1021     }
1022     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);
1023     if (X_t) {ierr = PetscMalloc1(numFaces*totDimBd,&u_t);CHKERRQ(ierr);}
1024     for (p = 0, f = 0; p < numPoints; ++p) {
1025       const PetscInt point = points[p];
1026       PetscScalar   *x     = NULL;
1027       PetscInt       i;
1028 
1029       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1030       if (dep != dim-1) continue;
1031       ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
1032       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
1033       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
1034       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1035       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1036       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1037       if (X_t) {
1038         ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1039         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1040         ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1041       }
1042       ++f;
1043     }
1044     for (f = 0; f < Nf; ++f) {
1045       PetscFE  fe;
1046       PetscInt numQuadPoints, Nb;
1047       /* Conforming batches */
1048       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1049       /* Remainder */
1050       PetscInt Nr, offset;
1051 
1052       ierr = PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1053       ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1054       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1055       ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1056       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1057       blockSize = Nb*numQuadPoints;
1058       batchSize = numBlocks * blockSize;
1059       ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1060       numChunks = numFaces / (numBatches*batchSize);
1061       Ne        = numChunks*numBatches*batchSize;
1062       Nr        = numFaces % (numBatches*batchSize);
1063       offset    = numFaces - Nr;
1064       geom.v0   = v0;
1065       geom.n    = n;
1066       geom.J    = J;
1067       geom.invJ = invJ;
1068       geom.detJ = detJ;
1069       ierr = PetscFEIntegrateBdResidual(fe, prob, f, Ne, geom, u, u_t, NULL, NULL, elemVec);CHKERRQ(ierr);
1070       geom.v0   = &v0[offset*dim];
1071       geom.n    = &n[offset*dim];
1072       geom.J    = &J[offset*dim*dim];
1073       geom.invJ = &invJ[offset*dim*dim];
1074       geom.detJ = &detJ[offset];
1075       ierr = PetscFEIntegrateBdResidual(fe, prob, f, Nr, geom, &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);CHKERRQ(ierr);
1076     }
1077     for (p = 0, f = 0; p < numPoints; ++p) {
1078       const PetscInt point = points[p];
1079 
1080       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1081       if (dep != dim-1) continue;
1082       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);CHKERRQ(ierr);}
1083       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*totDimBd], ADD_VALUES);CHKERRQ(ierr);
1084       ++f;
1085     }
1086     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1087     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1088     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1089     if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);}
1090   }
1091   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
1092   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
1093   PetscFunctionReturn(0);
1094 }
1095 
1096 #undef __FUNCT__
1097 #define __FUNCT__ "DMPlexComputeResidualFEM_Check_Internal"
1098 static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
1099 {
1100   DM                dmCh, dmAux;
1101   Vec               A;
1102   PetscDS           prob, probCh, probAux = NULL;
1103   PetscQuadrature   q;
1104   PetscCellGeometry geom;
1105   PetscSection      section, sectionAux;
1106   PetscReal        *v0, *J, *invJ, *detJ;
1107   PetscScalar      *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1108   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
1109   PetscInt          totDim, totDimAux, diffCell = 0;
1110   PetscErrorCode    ierr;
1111 
1112   PetscFunctionBegin;
1113   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1114   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1115   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1116   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1117   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1118   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1119   numCells = cEnd - cStart;
1120   ierr = PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);CHKERRQ(ierr);
1121   ierr = DMGetDS(dmCh, &probCh);CHKERRQ(ierr);
1122   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1123   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1124   if (dmAux) {
1125     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1126     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1127     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1128   }
1129   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
1130   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
1131   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);
1132   ierr = PetscMalloc1(numCells*totDim,&elemVecCh);CHKERRQ(ierr);
1133   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1134   for (c = cStart; c < cEnd; ++c) {
1135     PetscScalar *x = NULL, *x_t = NULL;
1136     PetscInt     i;
1137 
1138     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1139     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1140     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1141     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1142     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1143     if (X_t) {
1144       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1145       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1146       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1147     }
1148     if (dmAux) {
1149       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1150       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1151       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1152     }
1153   }
1154   for (f = 0; f < Nf; ++f) {
1155     PetscFE  fe, feCh;
1156     PetscInt numQuadPoints, Nb;
1157     /* Conforming batches */
1158     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1159     /* Remainder */
1160     PetscInt Nr, offset;
1161 
1162     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1163     ierr = PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);CHKERRQ(ierr);
1164     ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1165     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1166     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1167     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1168     blockSize = Nb*numQuadPoints;
1169     batchSize = numBlocks * blockSize;
1170     ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1171     numChunks = numCells / (numBatches*batchSize);
1172     Ne        = numChunks*numBatches*batchSize;
1173     Nr        = numCells % (numBatches*batchSize);
1174     offset    = numCells - Nr;
1175     geom.v0   = v0;
1176     geom.J    = J;
1177     geom.invJ = invJ;
1178     geom.detJ = detJ;
1179     ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, geom, u, u_t, probAux, a, elemVec);CHKERRQ(ierr);
1180     ierr = PetscFEIntegrateResidual(feCh, prob, f, Ne, geom, u, u_t, probAux, a, elemVecCh);CHKERRQ(ierr);
1181     geom.v0   = &v0[offset*dim];
1182     geom.J    = &J[offset*dim*dim];
1183     geom.invJ = &invJ[offset*dim*dim];
1184     geom.detJ = &detJ[offset];
1185     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);
1186     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);
1187   }
1188   for (c = cStart; c < cEnd; ++c) {
1189     PetscBool diff = PETSC_FALSE;
1190     PetscInt  d;
1191 
1192     for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1193     if (diff) {
1194       ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);CHKERRQ(ierr);
1195       ierr = DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);CHKERRQ(ierr);
1196       ierr = DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);CHKERRQ(ierr);
1197       ++diffCell;
1198     }
1199     if (diffCell > 9) break;
1200     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);CHKERRQ(ierr);
1201   }
1202   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1203   ierr = PetscFree(elemVecCh);CHKERRQ(ierr);
1204   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1205   PetscFunctionReturn(0);
1206 }
1207 
1208 #undef __FUNCT__
1209 #define __FUNCT__ "DMPlexSNESComputeResidualFEM"
1210 /*@
1211   DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
1212 
1213   Input Parameters:
1214 + dm - The mesh
1215 . X  - Local solution
1216 - user - The user context
1217 
1218   Output Parameter:
1219 . F  - Local output vector
1220 
1221   Level: developer
1222 
1223 .seealso: DMPlexComputeJacobianActionFEM()
1224 @*/
1225 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1226 {
1227   PetscObject    check;
1228   PetscErrorCode ierr;
1229 
1230   PetscFunctionBegin;
1231   /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
1232   ierr = PetscObjectQuery((PetscObject) dm, "dmCh", &check);CHKERRQ(ierr);
1233   if (check) {ierr = DMPlexComputeResidualFEM_Check_Internal(dm, X, NULL, F, user);CHKERRQ(ierr);}
1234   else       {ierr = DMPlexComputeResidualFEM_Internal(dm, X, NULL, F, user);CHKERRQ(ierr);}
1235   PetscFunctionReturn(0);
1236 }
1237 
1238 #undef __FUNCT__
1239 #define __FUNCT__ "DMPlexComputeJacobianFEM_Internal"
1240 PetscErrorCode DMPlexComputeJacobianFEM_Internal(DM dm, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user)
1241 {
1242   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1243   const char       *name  = "Jacobian";
1244   DM                dmAux;
1245   DMLabel           depth;
1246   Vec               A;
1247   PetscDS           prob, probAux = NULL;
1248   PetscQuadrature   quad;
1249   PetscCellGeometry geom;
1250   PetscSection      section, globalSection, sectionAux;
1251   PetscReal        *v0, *J, *invJ, *detJ;
1252   PetscScalar      *elemMat, *u, *u_t, *a = NULL;
1253   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1254   PetscInt          totDim, totDimBd, totDimAux, numBd, bd;
1255   PetscBool         isShell;
1256   PetscErrorCode    ierr;
1257 
1258   PetscFunctionBegin;
1259   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1260   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1261   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1262   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1263   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1264   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1265   ierr = PetscDSGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr);
1266   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1267   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1268   numCells = cEnd - cStart;
1269   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1270   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1271   if (dmAux) {
1272     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1273     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1274     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1275   }
1276   ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr);
1277   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1278   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);
1279   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1280   for (c = cStart; c < cEnd; ++c) {
1281     PetscScalar *x = NULL,  *x_t = NULL;
1282     PetscInt     i;
1283 
1284     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1285     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1286     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1287     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1288     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1289     if (X_t) {
1290       ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1291       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1292       ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
1293     }
1294     if (dmAux) {
1295       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1296       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1297       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1298     }
1299   }
1300   ierr = PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));CHKERRQ(ierr);
1301   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1302     PetscFE  fe;
1303     PetscInt numQuadPoints, Nb;
1304     /* Conforming batches */
1305     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1306     /* Remainder */
1307     PetscInt Nr, offset;
1308 
1309     ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1310     ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1311     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1312     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1313     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1314     blockSize = Nb*numQuadPoints;
1315     batchSize = numBlocks * blockSize;
1316     ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1317     numChunks = numCells / (numBatches*batchSize);
1318     Ne        = numChunks*numBatches*batchSize;
1319     Nr        = numCells % (numBatches*batchSize);
1320     offset    = numCells - Nr;
1321     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1322       geom.v0   = v0;
1323       geom.J    = J;
1324       geom.invJ = invJ;
1325       geom.detJ = detJ;
1326       ierr = PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, probAux, a, elemMat);CHKERRQ(ierr);
1327       geom.v0   = &v0[offset*dim];
1328       geom.J    = &J[offset*dim*dim];
1329       geom.invJ = &invJ[offset*dim*dim];
1330       geom.detJ = &detJ[offset];
1331       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);
1332     }
1333   }
1334   for (c = cStart; c < cEnd; ++c) {
1335     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[c*totDim*totDim]);CHKERRQ(ierr);}
1336     ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*totDim*totDim], ADD_VALUES);CHKERRQ(ierr);
1337   }
1338   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1339   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1340   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
1341   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
1342   ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
1343   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
1344   for (bd = 0; bd < numBd; ++bd) {
1345     const char     *bdLabel;
1346     DMLabel         label;
1347     IS              pointIS;
1348     const PetscInt *points;
1349     const PetscInt *values;
1350     PetscReal      *n;
1351     PetscInt        field, numValues, numPoints, p, dep, numFaces;
1352     PetscBool       isEssential;
1353 
1354     ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr);
1355     if (isEssential) continue;
1356     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1357     ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);
1358     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
1359     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
1360     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
1361     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1362       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1363       if (dep == dim-1) ++numFaces;
1364     }
1365     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);
1366     if (X_t) {ierr = PetscMalloc1(numFaces*totDimBd,&u_t);CHKERRQ(ierr);}
1367     for (p = 0, f = 0; p < numPoints; ++p) {
1368       const PetscInt point = points[p];
1369       PetscScalar   *x     = NULL;
1370       PetscInt       i;
1371 
1372       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
1373       if (dep != dim-1) continue;
1374       ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
1375       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
1376       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
1377       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1378       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1379       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
1380       if (X_t) {
1381         ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1382         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1383         ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
1384       }
1385       ++f;
1386     }
1387     ierr = PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));CHKERRQ(ierr);
1388     for (fieldI = 0; fieldI < Nf; ++fieldI) {
1389       PetscFE  fe;
1390       PetscInt numQuadPoints, Nb;
1391       /* Conforming batches */
1392       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1393       /* Remainder */
1394       PetscInt Nr, offset;
1395 
1396       ierr = PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1397       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1398       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1399       ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1400       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1401       blockSize = Nb*numQuadPoints;
1402       batchSize = numBlocks * blockSize;
1403       ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1404       numChunks = numFaces / (numBatches*batchSize);
1405       Ne        = numChunks*numBatches*batchSize;
1406       Nr        = numFaces % (numBatches*batchSize);
1407       offset    = numFaces - Nr;
1408       for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1409         geom.v0   = v0;
1410         geom.n    = n;
1411         geom.J    = J;
1412         geom.invJ = invJ;
1413         geom.detJ = detJ;
1414         ierr = PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, NULL, NULL, elemMat);CHKERRQ(ierr);
1415         geom.v0   = &v0[offset*dim];
1416         geom.n    = &n[offset*dim];
1417         geom.J    = &J[offset*dim*dim];
1418         geom.invJ = &invJ[offset*dim*dim];
1419         geom.detJ = &detJ[offset];
1420         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);
1421       }
1422     }
1423     for (p = 0, f = 0; p < numPoints; ++p) {
1424       const PetscInt point = points[p];
1425 
1426       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1427       if (dep != dim-1) continue;
1428       if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);CHKERRQ(ierr);}
1429       ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);CHKERRQ(ierr);
1430       ++f;
1431     }
1432     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1433     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1434     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1435     if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);}
1436   }
1437   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1438   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1439   if (mesh->printFEM) {
1440     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1441     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
1442     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1443   }
1444   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1445   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
1446   if (isShell) {
1447     JacActionCtx *jctx;
1448 
1449     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1450     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
1451   }
1452   PetscFunctionReturn(0);
1453 }
1454 
1455 #undef __FUNCT__
1456 #define __FUNCT__ "DMPlexSNESComputeJacobianFEM"
1457 /*@
1458   DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1459 
1460   Input Parameters:
1461 + dm - The mesh
1462 . X  - Local input vector
1463 - user - The user context
1464 
1465   Output Parameter:
1466 . Jac  - Jacobian matrix
1467 
1468   Note:
1469   The first member of the user context must be an FEMContext.
1470 
1471   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1472   like a GPU, or vectorize on a multicore machine.
1473 
1474   Level: developer
1475 
1476 .seealso: FormFunctionLocal()
1477 @*/
1478 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1479 {
1480   PetscErrorCode ierr;
1481 
1482   PetscFunctionBegin;
1483   ierr = DMPlexComputeJacobianFEM_Internal(dm, X, NULL, Jac, JacP, user);CHKERRQ(ierr);
1484   PetscFunctionReturn(0);
1485 }
1486