xref: /petsc/src/snes/utils/dmplexsnes.c (revision 8ad47952706b5b3a4d2ddd55b418d583f25bc7ec) 
1 #include <petscdmplex.h> /*I "petscdmplex.h" I*/
2 #include <petscsnes.h>      /*I "petscsnes.h" I*/
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "DMInterpolationCreate"
6 PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx) {
7   PetscErrorCode ierr;
8 
9   PetscFunctionBegin;
10   PetscValidPointer(ctx, 2);
11   ierr = PetscMalloc(sizeof(struct _DMInterpolationInfo), ctx);CHKERRQ(ierr);
12   (*ctx)->comm   = comm;
13   (*ctx)->dim    = -1;
14   (*ctx)->nInput = 0;
15   (*ctx)->points = PETSC_NULL;
16   (*ctx)->cells  = PETSC_NULL;
17   (*ctx)->n      = -1;
18   (*ctx)->coords = PETSC_NULL;
19   PetscFunctionReturn(0);
20 }
21 
22 #undef __FUNCT__
23 #define __FUNCT__ "DMInterpolationSetDim"
24 PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim) {
25   PetscFunctionBegin;
26   if ((dim < 1) || (dim > 3)) {SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);}
27   ctx->dim = dim;
28   PetscFunctionReturn(0);
29 }
30 
31 #undef __FUNCT__
32 #define __FUNCT__ "DMInterpolationGetDim"
33 PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim) {
34   PetscFunctionBegin;
35   PetscValidIntPointer(dim, 2);
36   *dim = ctx->dim;
37   PetscFunctionReturn(0);
38 }
39 
40 #undef __FUNCT__
41 #define __FUNCT__ "DMInterpolationSetDof"
42 PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof) {
43   PetscFunctionBegin;
44   if (dof < 1) {SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);}
45   ctx->dof = dof;
46   PetscFunctionReturn(0);
47 }
48 
49 #undef __FUNCT__
50 #define __FUNCT__ "DMInterpolationGetDof"
51 PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof) {
52   PetscFunctionBegin;
53   PetscValidIntPointer(dof, 2);
54   *dof = ctx->dof;
55   PetscFunctionReturn(0);
56 }
57 
58 #undef __FUNCT__
59 #define __FUNCT__ "DMInterpolationAddPoints"
60 PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[]) {
61   PetscErrorCode ierr;
62 
63   PetscFunctionBegin;
64   if (ctx->dim < 0) {
65     SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
66   }
67   if (ctx->points) {
68     SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
69   }
70   ctx->nInput = n;
71   ierr = PetscMalloc(n*ctx->dim * sizeof(PetscReal), &ctx->points);CHKERRQ(ierr);
72   ierr = PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));CHKERRQ(ierr);
73   PetscFunctionReturn(0);
74 }
75 
76 #undef __FUNCT__
77 #define __FUNCT__ "DMInterpolationSetUp"
78 PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints) {
79   MPI_Comm       comm = ctx->comm;
80   PetscScalar   *a;
81   PetscInt       p, q, i;
82   PetscMPIInt    rank, size;
83   PetscErrorCode ierr;
84 
85   PetscFunctionBegin;
86   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
87   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
88   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
89   if (ctx->dim < 0) {
90     SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
91   }
92   // Locate points
93   Vec             pointVec;
94   IS              cellIS;
95   PetscLayout     layout;
96   PetscReal      *globalPoints;
97   const PetscInt *ranges;
98   PetscMPIInt    *counts, *displs;
99   const PetscInt *foundCells;
100   PetscMPIInt    *foundProcs, *globalProcs;
101   PetscInt        n = ctx->nInput, N;
102 
103   if (!redundantPoints) {
104     ierr = PetscLayoutCreate(comm, &layout);CHKERRQ(ierr);
105     ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr);
106     ierr = PetscLayoutSetLocalSize(layout, n);CHKERRQ(ierr);
107     ierr = PetscLayoutSetUp(layout);CHKERRQ(ierr);
108     ierr = PetscLayoutGetSize(layout, &N);CHKERRQ(ierr);
109     /* Communicate all points to all processes */
110     ierr = PetscMalloc3(N*ctx->dim,PetscReal,&globalPoints,size,PetscMPIInt,&counts,size,PetscMPIInt,&displs);CHKERRQ(ierr);
111     ierr = PetscLayoutGetRanges(layout, &ranges);CHKERRQ(ierr);
112     for (p = 0; p < size; ++p) {
113       counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
114       displs[p] = ranges[p]*ctx->dim;
115     }
116     ierr = MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);CHKERRQ(ierr);
117   } else {
118     N = n;
119     globalPoints = ctx->points;
120   }
121 #if 0
122   ierr = PetscMalloc3(N,PetscInt,&foundCells,N,PetscMPIInt,&foundProcs,N,PetscMPIInt,&globalProcs);CHKERRQ(ierr);
123   //foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]);
124 #else
125   ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N, globalPoints, &pointVec);CHKERRQ(ierr);
126   ierr = PetscMalloc2(N,PetscMPIInt,&foundProcs,N,PetscMPIInt,&globalProcs);CHKERRQ(ierr);
127   ierr = DMLocatePoints(dm, pointVec, &cellIS);CHKERRQ(ierr);
128   ierr = ISGetIndices(cellIS, &foundCells);CHKERRQ(ierr);
129 #endif
130   for (p = 0; p < N; ++p) {
131     if (foundCells[p] >= 0) {
132       foundProcs[p] = rank;
133     } else {
134       foundProcs[p] = size;
135     }
136   }
137   /* Let the lowest rank process own each point */
138   ierr = MPI_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);CHKERRQ(ierr);
139   ctx->n = 0;
140   for (p = 0; p < N; ++p) {
141     if (globalProcs[p] == size) {
142       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);
143     } else if (globalProcs[p] == rank) {
144       ctx->n++;
145     }
146   }
147   /* Create coordinates vector and array of owned cells */
148   ierr = PetscMalloc(ctx->n * sizeof(PetscInt), &ctx->cells);CHKERRQ(ierr);
149   ierr = VecCreate(comm, &ctx->coords);CHKERRQ(ierr);
150   ierr = VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);CHKERRQ(ierr);
151   ierr = VecSetBlockSize(ctx->coords, ctx->dim);CHKERRQ(ierr);
152   ierr = VecSetFromOptions(ctx->coords);CHKERRQ(ierr);
153   ierr = VecGetArray(ctx->coords, &a);CHKERRQ(ierr);
154   for (p = 0, q = 0, i = 0; p < N; ++p) {
155     if (globalProcs[p] == rank) {
156       PetscInt d;
157 
158       for (d = 0; d < ctx->dim; ++d, ++i) {
159         a[i] = globalPoints[p*ctx->dim+d];
160       }
161       ctx->cells[q++] = foundCells[p];
162     }
163   }
164   ierr = VecRestoreArray(ctx->coords, &a);CHKERRQ(ierr);
165 #if 0
166   ierr = PetscFree3(foundCells,foundProcs,globalProcs);CHKERRQ(ierr);
167 #else
168   ierr = PetscFree2(foundProcs,globalProcs);CHKERRQ(ierr);
169   ierr = ISRestoreIndices(cellIS, &foundCells);CHKERRQ(ierr);
170   ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
171   ierr = VecDestroy(&pointVec);CHKERRQ(ierr);
172 #endif
173   if (!redundantPoints) {
174     ierr = PetscFree3(globalPoints,counts,displs);CHKERRQ(ierr);
175     ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr);
176   }
177   PetscFunctionReturn(0);
178 }
179 
180 #undef __FUNCT__
181 #define __FUNCT__ "DMInterpolationGetCoordinates"
182 PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates) {
183   PetscFunctionBegin;
184   PetscValidPointer(coordinates, 2);
185   if (!ctx->coords) {SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");}
186   *coordinates = ctx->coords;
187   PetscFunctionReturn(0);
188 }
189 
190 #undef __FUNCT__
191 #define __FUNCT__ "DMInterpolationGetVector"
192 PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v) {
193   PetscErrorCode ierr;
194 
195   PetscFunctionBegin;
196   PetscValidPointer(v, 2);
197   if (!ctx->coords) {SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");}
198   ierr = VecCreate(ctx->comm, v);CHKERRQ(ierr);
199   ierr = VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);CHKERRQ(ierr);
200   ierr = VecSetBlockSize(*v, ctx->dof);CHKERRQ(ierr);
201   ierr = VecSetFromOptions(*v);CHKERRQ(ierr);
202   PetscFunctionReturn(0);
203 }
204 
205 #undef __FUNCT__
206 #define __FUNCT__ "DMInterpolationRestoreVector"
207 PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v) {
208   PetscErrorCode ierr;
209 
210   PetscFunctionBegin;
211   PetscValidPointer(v, 2);
212   if (!ctx->coords) {SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");}
213   ierr = VecDestroy(v);CHKERRQ(ierr);
214   PetscFunctionReturn(0);
215 }
216 
217 #undef __FUNCT__
218 #define __FUNCT__ "DMInterpolate_Simplex_Private"
219 PetscErrorCode DMInterpolate_Simplex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v) {
220   PetscReal     *v0, *J, *invJ, detJ;
221   PetscScalar   *a, *coords;
222   PetscInt       p;
223   PetscErrorCode ierr;
224 
225   PetscFunctionBegin;
226   ierr = PetscMalloc3(ctx->dim,PetscReal,&v0,ctx->dim*ctx->dim,PetscReal,&J,ctx->dim*ctx->dim,PetscReal,&invJ);CHKERRQ(ierr);
227   ierr = VecGetArray(ctx->coords, &coords);CHKERRQ(ierr);
228   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
229   for(p = 0; p < ctx->n; ++p) {
230     PetscInt           c = ctx->cells[p];
231     const PetscScalar *x;
232     PetscReal          xi[4];
233     PetscInt           d, f, comp;
234 
235     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
236     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
237     ierr = DMPlexVecGetClosure(dm, PETSC_NULL, xLocal, c, PETSC_NULL, &x);CHKERRQ(ierr);
238     for(comp = 0; comp < ctx->dof; ++comp) {
239       a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
240     }
241     for(d = 0; d < ctx->dim; ++d) {
242       xi[d] = 0.0;
243       for(f = 0; f < ctx->dim; ++f) {
244         xi[d] += invJ[d*ctx->dim+f]*0.5*(coords[p*ctx->dim+f] - v0[f]);
245       }
246       for(comp = 0; comp < ctx->dof; ++comp) {
247         a[p*ctx->dof+comp] += (x[(d+1)*ctx->dof+comp] - x[0*ctx->dof+comp])*xi[d];
248       }
249     }
250     ierr = DMPlexVecRestoreClosure(dm, PETSC_NULL, xLocal, c, PETSC_NULL, &x);CHKERRQ(ierr);
251   }
252   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
253   ierr = VecRestoreArray(ctx->coords, &coords);CHKERRQ(ierr);
254   ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr);
255   PetscFunctionReturn(0);
256 }
257 
258 #undef __FUNCT__
259 #define __FUNCT__ "QuadMap_Private"
260 PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
261 {
262   const PetscScalar*vertices = (const PetscScalar *) ctx;
263   const PetscScalar x0   = vertices[0];
264   const PetscScalar y0   = vertices[1];
265   const PetscScalar x1   = vertices[2];
266   const PetscScalar y1   = vertices[3];
267   const PetscScalar x2   = vertices[4];
268   const PetscScalar y2   = vertices[5];
269   const PetscScalar x3   = vertices[6];
270   const PetscScalar y3   = vertices[7];
271   const PetscScalar f_1  = x1 - x0;
272   const PetscScalar g_1  = y1 - y0;
273   const PetscScalar f_3  = x3 - x0;
274   const PetscScalar g_3  = y3 - y0;
275   const PetscScalar f_01 = x2 - x1 - x3 + x0;
276   const PetscScalar g_01 = y2 - y1 - y3 + y0;
277   PetscScalar      *ref, *real;
278   PetscErrorCode    ierr;
279 
280   PetscFunctionBegin;
281   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
282   ierr = VecGetArray(Xreal, &real);CHKERRQ(ierr);
283   {
284     const PetscScalar p0 = ref[0];
285     const PetscScalar p1 = ref[1];
286 
287     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
288     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
289   }
290   ierr = PetscLogFlops(28);CHKERRQ(ierr);
291   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
292   ierr = VecRestoreArray(Xreal, &real);CHKERRQ(ierr);
293   PetscFunctionReturn(0);
294 }
295 
296 #undef __FUNCT__
297 #define __FUNCT__ "QuadJacobian_Private"
298 PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat *J, Mat *M, MatStructure *flag, void *ctx)
299 {
300   const PetscScalar*vertices = (const PetscScalar *) ctx;
301   const PetscScalar x0   = vertices[0];
302   const PetscScalar y0   = vertices[1];
303   const PetscScalar x1   = vertices[2];
304   const PetscScalar y1   = vertices[3];
305   const PetscScalar x2   = vertices[4];
306   const PetscScalar y2   = vertices[5];
307   const PetscScalar x3   = vertices[6];
308   const PetscScalar y3   = vertices[7];
309   const PetscScalar f_01 = x2 - x1 - x3 + x0;
310   const PetscScalar g_01 = y2 - y1 - y3 + y0;
311   PetscScalar      *ref;
312   PetscErrorCode    ierr;
313 
314   PetscFunctionBegin;
315   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
316   {
317     const PetscScalar x = ref[0];
318     const PetscScalar y = ref[1];
319     const PetscInt    rows[2]   = {0, 1};
320     const PetscScalar values[4] = {(x1 - x0 + f_01*y) * 0.5, (x3 - x0 + f_01*x) * 0.5,
321                                    (y1 - y0 + g_01*y) * 0.5, (y3 - y0 + g_01*x) * 0.5};
322     ierr = MatSetValues(*J, 2, rows, 2, rows, values, INSERT_VALUES);CHKERRQ(ierr);
323   }
324   ierr = PetscLogFlops(30);CHKERRQ(ierr);
325   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
326   ierr = MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
327   ierr = MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
328   PetscFunctionReturn(0);
329 }
330 
331 #undef __FUNCT__
332 #define __FUNCT__ "DMInterpolate_Quad_Private"
333 PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v) {
334   SNES           snes;
335   KSP            ksp;
336   PC             pc;
337   Vec            coordsLocal, r, ref, real;
338   Mat            J;
339   PetscScalar   *a, *coords;
340   PetscInt       p;
341   PetscErrorCode ierr;
342 
343   PetscFunctionBegin;
344   ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
345   ierr = SNESCreate(PETSC_COMM_SELF, &snes);CHKERRQ(ierr);
346   ierr = SNESSetOptionsPrefix(snes, "quad_interp_");CHKERRQ(ierr);
347   ierr = VecCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
348   ierr = VecSetSizes(r, 2, 2);CHKERRQ(ierr);
349   ierr = VecSetFromOptions(r);CHKERRQ(ierr);
350   ierr = VecDuplicate(r, &ref);CHKERRQ(ierr);
351   ierr = VecDuplicate(r, &real);CHKERRQ(ierr);
352   ierr = MatCreate(PETSC_COMM_SELF, &J);CHKERRQ(ierr);
353   ierr = MatSetSizes(J, 2, 2, 2, 2);CHKERRQ(ierr);
354   ierr = MatSetType(J, MATSEQDENSE);CHKERRQ(ierr);
355   ierr = MatSetUp(J);CHKERRQ(ierr);
356   ierr = SNESSetFunction(snes, r, QuadMap_Private, PETSC_NULL);CHKERRQ(ierr);
357   ierr = SNESSetJacobian(snes, J, J, QuadJacobian_Private, PETSC_NULL);CHKERRQ(ierr);
358   ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
359   ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
360   ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
361   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
362 
363   ierr = VecGetArray(ctx->coords, &coords);CHKERRQ(ierr);
364   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
365   for (p = 0; p < ctx->n; ++p) {
366     const PetscScalar *x, *vertices;
367     PetscScalar       *xi;
368     PetscInt           c = ctx->cells[p], comp, coordSize, xSize;
369 
370     /* Can make this do all points at once */
371     ierr = DMPlexVecGetClosure(dm, PETSC_NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
372     if (4*2 != coordSize) {SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);}
373     ierr = DMPlexVecGetClosure(dm, PETSC_NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
374     if (4*ctx->dof != xSize) {SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 4*ctx->dof);}
375     ierr = SNESSetFunction(snes, PETSC_NULL, PETSC_NULL, (void *) vertices);CHKERRQ(ierr);
376     ierr = SNESSetJacobian(snes, PETSC_NULL, PETSC_NULL, PETSC_NULL, (void *) vertices);CHKERRQ(ierr);
377     ierr = VecGetArray(real, &xi);CHKERRQ(ierr);
378     xi[0] = coords[p*ctx->dim+0];
379     xi[1] = coords[p*ctx->dim+1];
380     ierr = VecRestoreArray(real, &xi);CHKERRQ(ierr);
381     ierr = SNESSolve(snes, real, ref);CHKERRQ(ierr);
382     ierr = VecGetArray(ref, &xi);CHKERRQ(ierr);
383     for (comp = 0; comp < ctx->dof; ++comp) {
384       a[p*ctx->dof+comp] = x[0*ctx->dof+comp]*(1 - xi[0])*(1 - xi[1]) + x[1*ctx->dof+comp]*xi[0]*(1 - xi[1]) + x[2*ctx->dof+comp]*xi[0]*xi[1] + x[3*ctx->dof+comp]*(1 - xi[0])*xi[1];
385     }
386     ierr = VecRestoreArray(ref, &xi);CHKERRQ(ierr);
387     ierr = DMPlexVecRestoreClosure(dm, PETSC_NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
388   }
389   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
390   ierr = VecRestoreArray(ctx->coords, &coords);CHKERRQ(ierr);
391 
392   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
393   ierr = VecDestroy(&r);CHKERRQ(ierr);
394   ierr = VecDestroy(&ref);CHKERRQ(ierr);
395   ierr = VecDestroy(&real);CHKERRQ(ierr);
396   ierr = MatDestroy(&J);CHKERRQ(ierr);
397   PetscFunctionReturn(0);
398 }
399 
400 #undef __FUNCT__
401 #define __FUNCT__ "HexMap_Private"
402 PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
403 {
404   const PetscScalar*vertices = (const PetscScalar *) ctx;
405   const PetscScalar x0 = vertices[0];
406   const PetscScalar y0 = vertices[1];
407   const PetscScalar z0 = vertices[2];
408   const PetscScalar x1 = vertices[3];
409   const PetscScalar y1 = vertices[4];
410   const PetscScalar z1 = vertices[5];
411   const PetscScalar x2 = vertices[6];
412   const PetscScalar y2 = vertices[7];
413   const PetscScalar z2 = vertices[8];
414   const PetscScalar x3 = vertices[9];
415   const PetscScalar y3 = vertices[10];
416   const PetscScalar z3 = vertices[11];
417   const PetscScalar x4 = vertices[12];
418   const PetscScalar y4 = vertices[13];
419   const PetscScalar z4 = vertices[14];
420   const PetscScalar x5 = vertices[15];
421   const PetscScalar y5 = vertices[16];
422   const PetscScalar z5 = vertices[17];
423   const PetscScalar x6 = vertices[18];
424   const PetscScalar y6 = vertices[19];
425   const PetscScalar z6 = vertices[20];
426   const PetscScalar x7 = vertices[21];
427   const PetscScalar y7 = vertices[22];
428   const PetscScalar z7 = vertices[23];
429   const PetscScalar f_1 = x1 - x0;
430   const PetscScalar g_1 = y1 - y0;
431   const PetscScalar h_1 = z1 - z0;
432   const PetscScalar f_3 = x3 - x0;
433   const PetscScalar g_3 = y3 - y0;
434   const PetscScalar h_3 = z3 - z0;
435   const PetscScalar f_4 = x4 - x0;
436   const PetscScalar g_4 = y4 - y0;
437   const PetscScalar h_4 = z4 - z0;
438   const PetscScalar f_01 = x2 - x1 - x3 + x0;
439   const PetscScalar g_01 = y2 - y1 - y3 + y0;
440   const PetscScalar h_01 = z2 - z1 - z3 + z0;
441   const PetscScalar f_12 = x7 - x3 - x4 + x0;
442   const PetscScalar g_12 = y7 - y3 - y4 + y0;
443   const PetscScalar h_12 = z7 - z3 - z4 + z0;
444   const PetscScalar f_02 = x5 - x1 - x4 + x0;
445   const PetscScalar g_02 = y5 - y1 - y4 + y0;
446   const PetscScalar h_02 = z5 - z1 - z4 + z0;
447   const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
448   const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
449   const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
450   PetscScalar      *ref, *real;
451   PetscErrorCode    ierr;
452 
453   PetscFunctionBegin;
454   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
455   ierr = VecGetArray(Xreal, &real);CHKERRQ(ierr);
456   {
457     const PetscScalar p0 = ref[0];
458     const PetscScalar p1 = ref[1];
459     const PetscScalar p2 = ref[2];
460 
461     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;
462     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;
463     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;
464   }
465   ierr = PetscLogFlops(114);CHKERRQ(ierr);
466   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
467   ierr = VecRestoreArray(Xreal, &real);CHKERRQ(ierr);
468   PetscFunctionReturn(0);
469 }
470 
471 #undef __FUNCT__
472 #define __FUNCT__ "HexJacobian_Private"
473 PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat *J, Mat *M, MatStructure *flag, void *ctx)
474 {
475   const PetscScalar*vertices = (const PetscScalar *) ctx;
476   const PetscScalar x0 = vertices[0];
477   const PetscScalar y0 = vertices[1];
478   const PetscScalar z0 = vertices[2];
479   const PetscScalar x1 = vertices[3];
480   const PetscScalar y1 = vertices[4];
481   const PetscScalar z1 = vertices[5];
482   const PetscScalar x2 = vertices[6];
483   const PetscScalar y2 = vertices[7];
484   const PetscScalar z2 = vertices[8];
485   const PetscScalar x3 = vertices[9];
486   const PetscScalar y3 = vertices[10];
487   const PetscScalar z3 = vertices[11];
488   const PetscScalar x4 = vertices[12];
489   const PetscScalar y4 = vertices[13];
490   const PetscScalar z4 = vertices[14];
491   const PetscScalar x5 = vertices[15];
492   const PetscScalar y5 = vertices[16];
493   const PetscScalar z5 = vertices[17];
494   const PetscScalar x6 = vertices[18];
495   const PetscScalar y6 = vertices[19];
496   const PetscScalar z6 = vertices[20];
497   const PetscScalar x7 = vertices[21];
498   const PetscScalar y7 = vertices[22];
499   const PetscScalar z7 = vertices[23];
500   const PetscScalar f_xy = x2 - x1 - x3 + x0;
501   const PetscScalar g_xy = y2 - y1 - y3 + y0;
502   const PetscScalar h_xy = z2 - z1 - z3 + z0;
503   const PetscScalar f_yz = x7 - x3 - x4 + x0;
504   const PetscScalar g_yz = y7 - y3 - y4 + y0;
505   const PetscScalar h_yz = z7 - z3 - z4 + z0;
506   const PetscScalar f_xz = x5 - x1 - x4 + x0;
507   const PetscScalar g_xz = y5 - y1 - y4 + y0;
508   const PetscScalar h_xz = z5 - z1 - z4 + z0;
509   const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
510   const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
511   const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
512   PetscScalar      *ref;
513   PetscErrorCode    ierr;
514 
515   PetscFunctionBegin;
516   ierr = VecGetArray(Xref,  &ref);CHKERRQ(ierr);
517   {
518     const PetscScalar x = ref[0];
519     const PetscScalar y = ref[1];
520     const PetscScalar z = ref[2];
521     const PetscInt    rows[3]   = {0, 1, 2};
522     const PetscScalar values[9] = {
523       (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0,
524       (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0,
525       (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0,
526       (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0,
527       (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0,
528       (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0,
529       (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0,
530       (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0,
531       (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0};
532     ierr = MatSetValues(*J, 3, rows, 3, rows, values, INSERT_VALUES);CHKERRQ(ierr);
533   }
534   ierr = PetscLogFlops(152);CHKERRQ(ierr);
535   ierr = VecRestoreArray(Xref,  &ref);CHKERRQ(ierr);
536   ierr = MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
537   ierr = MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
538   PetscFunctionReturn(0);
539 }
540 
541 #undef __FUNCT__
542 #define __FUNCT__ "DMInterpolate_Hex_Private"
543 PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v) {
544   SNES           snes;
545   KSP            ksp;
546   PC             pc;
547   Vec            coordsLocal, r, ref, real;
548   Mat            J;
549   PetscScalar   *a, *coords;
550   PetscInt       p;
551   PetscErrorCode ierr;
552 
553   PetscFunctionBegin;
554   ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
555   ierr = SNESCreate(PETSC_COMM_SELF, &snes);CHKERRQ(ierr);
556   ierr = SNESSetOptionsPrefix(snes, "hex_interp_");CHKERRQ(ierr);
557   ierr = VecCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
558   ierr = VecSetSizes(r, 3, 3);CHKERRQ(ierr);
559   ierr = VecSetFromOptions(r);CHKERRQ(ierr);
560   ierr = VecDuplicate(r, &ref);CHKERRQ(ierr);
561   ierr = VecDuplicate(r, &real);CHKERRQ(ierr);
562   ierr = MatCreate(PETSC_COMM_SELF, &J);CHKERRQ(ierr);
563   ierr = MatSetSizes(J, 3, 3, 3, 3);CHKERRQ(ierr);
564   ierr = MatSetType(J, MATSEQDENSE);CHKERRQ(ierr);
565   ierr = MatSetUp(J);CHKERRQ(ierr);
566   ierr = SNESSetFunction(snes, r, HexMap_Private, PETSC_NULL);CHKERRQ(ierr);
567   ierr = SNESSetJacobian(snes, J, J, HexJacobian_Private, PETSC_NULL);CHKERRQ(ierr);
568   ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
569   ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
570   ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
571   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
572 
573   ierr = VecGetArray(ctx->coords, &coords);CHKERRQ(ierr);
574   ierr = VecGetArray(v, &a);CHKERRQ(ierr);
575   for(p = 0; p < ctx->n; ++p) {
576     const PetscScalar *x, *vertices;
577     PetscScalar       *xi;
578     PetscInt           c = ctx->cells[p], comp, coordSize, xSize;
579 
580     /* Can make this do all points at once */
581     ierr = DMPlexVecGetClosure(dm, PETSC_NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
582     if (8*3 != coordSize) {SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);}
583     ierr = DMPlexVecGetClosure(dm, PETSC_NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
584     if (8*ctx->dof != xSize) {SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);}
585     ierr = SNESSetFunction(snes, PETSC_NULL, PETSC_NULL, (void *) vertices);CHKERRQ(ierr);
586     ierr = SNESSetJacobian(snes, PETSC_NULL, PETSC_NULL, PETSC_NULL, (void *) vertices);CHKERRQ(ierr);
587     ierr = VecGetArray(real, &xi);CHKERRQ(ierr);
588     xi[0] = coords[p*ctx->dim+0];
589     xi[1] = coords[p*ctx->dim+1];
590     xi[2] = coords[p*ctx->dim+2];
591     ierr = VecRestoreArray(real, &xi);CHKERRQ(ierr);
592     ierr = SNESSolve(snes, real, ref);CHKERRQ(ierr);
593     ierr = VecGetArray(ref, &xi);CHKERRQ(ierr);
594     for (comp = 0; comp < ctx->dof; ++comp) {
595       a[p*ctx->dof+comp] =
596         x[0*ctx->dof+comp]*(1-xi[0])*(1-xi[1])*(1-xi[2]) +
597         x[1*ctx->dof+comp]*    xi[0]*(1-xi[1])*(1-xi[2]) +
598         x[2*ctx->dof+comp]*    xi[0]*    xi[1]*(1-xi[2]) +
599         x[3*ctx->dof+comp]*(1-xi[0])*    xi[1]*(1-xi[2]) +
600         x[4*ctx->dof+comp]*(1-xi[0])*(1-xi[1])*   xi[2] +
601         x[5*ctx->dof+comp]*    xi[0]*(1-xi[1])*   xi[2] +
602         x[6*ctx->dof+comp]*    xi[0]*    xi[1]*   xi[2] +
603         x[7*ctx->dof+comp]*(1-xi[0])*    xi[1]*   xi[2];
604     }
605     ierr = VecRestoreArray(ref, &xi);CHKERRQ(ierr);
606     ierr = DMPlexVecRestoreClosure(dm, PETSC_NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
607     ierr = DMPlexVecRestoreClosure(dm, PETSC_NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
608   }
609   ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
610   ierr = VecRestoreArray(ctx->coords, &coords);CHKERRQ(ierr);
611 
612   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
613   ierr = VecDestroy(&r);CHKERRQ(ierr);
614   ierr = VecDestroy(&ref);CHKERRQ(ierr);
615   ierr = VecDestroy(&real);CHKERRQ(ierr);
616   ierr = MatDestroy(&J);CHKERRQ(ierr);
617   PetscFunctionReturn(0);
618 }
619 
620 #undef __FUNCT__
621 #define __FUNCT__ "DMInterpolationEvaluate"
622 /*
623   Input Parameters:
624 + ctx - The DMInterpolationInfo context
625 . dm  - The DM
626 - x   - The local vector containing the field to be interpolated
627 
628   Output Parameters:
629 . v   - The vector containing the interpolated values
630 */
631 PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v) {
632   PetscInt       dim, coneSize, n;
633   PetscErrorCode ierr;
634 
635   PetscFunctionBegin;
636   PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
637   PetscValidHeaderSpecific(x, VEC_CLASSID, 3);
638   PetscValidHeaderSpecific(v, VEC_CLASSID, 4);
639   ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
640   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);}
641   if (n) {
642     ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
643     ierr = DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);CHKERRQ(ierr);
644     if (dim == 2) {
645       if (coneSize == 3) {
646         ierr = DMInterpolate_Simplex_Private(ctx, dm, x, v);CHKERRQ(ierr);
647       } else if (coneSize == 4) {
648         ierr = DMInterpolate_Quad_Private(ctx, dm, x, v);CHKERRQ(ierr);
649       } else {
650         SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
651       }
652     } else if (dim == 3) {
653       if (coneSize == 4) {
654         ierr = DMInterpolate_Simplex_Private(ctx, dm, x, v);CHKERRQ(ierr);
655       } else {
656         ierr = DMInterpolate_Hex_Private(ctx, dm, x, v);CHKERRQ(ierr);
657       }
658     } else {
659       SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
660     }
661   }
662   PetscFunctionReturn(0);
663 }
664 
665 #undef __FUNCT__
666 #define __FUNCT__ "DMInterpolationDestroy"
667 PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx) {
668   PetscErrorCode ierr;
669 
670   PetscFunctionBegin;
671   PetscValidPointer(ctx, 2);
672   ierr = VecDestroy(&(*ctx)->coords);CHKERRQ(ierr);
673   ierr = PetscFree((*ctx)->points);CHKERRQ(ierr);
674   ierr = PetscFree((*ctx)->cells);CHKERRQ(ierr);
675   ierr = PetscFree(*ctx);CHKERRQ(ierr);
676   *ctx = PETSC_NULL;
677   PetscFunctionReturn(0);
678 }
679