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