Actual source code: dmplexsnes.c
petsc-3.5.4 2015-05-23
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: PetscMalloc1(n*ctx->dim, &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,&globalPoints,size,&counts,size,&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;
124: globalPoints = ctx->points;
125: counts = displs = NULL;
126: layout = NULL;
127: }
128: #if 0
129: PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
130: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
131: #else
132: #if defined(PETSC_USE_COMPLEX)
133: PetscMalloc1(N,&globalPointsScalar);
134: for (i=0; i<N; i++) globalPointsScalar[i] = globalPoints[i];
135: #else
136: globalPointsScalar = globalPoints;
137: #endif
138: VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
139: PetscMalloc2(N,&foundProcs,N,&globalProcs);
140: DMLocatePoints(dm, pointVec, &cellIS);
141: ISGetIndices(cellIS, &foundCells);
142: #endif
143: for (p = 0; p < N; ++p) {
144: if (foundCells[p] >= 0) foundProcs[p] = rank;
145: else foundProcs[p] = size;
146: }
147: /* Let the lowest rank process own each point */
148: MPI_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
149: ctx->n = 0;
150: for (p = 0; p < N; ++p) {
151: 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);
152: else if (globalProcs[p] == rank) ctx->n++;
153: }
154: /* Create coordinates vector and array of owned cells */
155: PetscMalloc1(ctx->n, &ctx->cells);
156: VecCreate(comm, &ctx->coords);
157: VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
158: VecSetBlockSize(ctx->coords, ctx->dim);
159: VecSetType(ctx->coords,VECSTANDARD);
160: VecGetArray(ctx->coords, &a);
161: for (p = 0, q = 0, i = 0; p < N; ++p) {
162: if (globalProcs[p] == rank) {
163: PetscInt d;
165: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
166: ctx->cells[q++] = foundCells[p];
167: }
168: }
169: VecRestoreArray(ctx->coords, &a);
170: #if 0
171: PetscFree3(foundCells,foundProcs,globalProcs);
172: #else
173: PetscFree2(foundProcs,globalProcs);
174: ISRestoreIndices(cellIS, &foundCells);
175: ISDestroy(&cellIS);
176: VecDestroy(&pointVec);
177: #endif
178: if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
179: if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
180: PetscLayoutDestroy(&layout);
181: return(0);
182: }
186: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
187: {
190: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
191: *coordinates = ctx->coords;
192: return(0);
193: }
197: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
198: {
203: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
204: VecCreate(ctx->comm, v);
205: VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
206: VecSetBlockSize(*v, ctx->dof);
207: VecSetType(*v,VECSTANDARD);
208: return(0);
209: }
213: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
214: {
219: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
220: VecDestroy(v);
221: return(0);
222: }
226: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
227: {
228: PetscReal *v0, *J, *invJ, detJ;
229: PetscScalar *a, *coords;
230: PetscInt p;
234: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
235: VecGetArray(ctx->coords, &coords);
236: VecGetArray(v, &a);
237: for (p = 0; p < ctx->n; ++p) {
238: PetscInt c = ctx->cells[p];
239: PetscScalar *x = NULL;
240: PetscReal xi[4];
241: PetscInt d, f, comp;
243: DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
244: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
245: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
246: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
248: for (d = 0; d < ctx->dim; ++d) {
249: xi[d] = 0.0;
250: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
251: 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];
252: }
253: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
254: }
255: VecRestoreArray(v, &a);
256: VecRestoreArray(ctx->coords, &coords);
257: PetscFree3(v0, J, invJ);
258: return(0);
259: }
263: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
264: {
265: PetscReal *v0, *J, *invJ, detJ;
266: PetscScalar *a, *coords;
267: PetscInt p;
271: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
272: VecGetArray(ctx->coords, &coords);
273: VecGetArray(v, &a);
274: for (p = 0; p < ctx->n; ++p) {
275: PetscInt c = ctx->cells[p];
276: const PetscInt order[3] = {2, 1, 3};
277: PetscScalar *x = NULL;
278: PetscReal xi[4];
279: PetscInt d, f, comp;
281: DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
282: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
283: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
284: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
286: for (d = 0; d < ctx->dim; ++d) {
287: xi[d] = 0.0;
288: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
289: 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];
290: }
291: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
292: }
293: VecRestoreArray(v, &a);
294: VecRestoreArray(ctx->coords, &coords);
295: PetscFree3(v0, J, invJ);
296: return(0);
297: }
301: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
302: {
303: const PetscScalar *vertices = (const PetscScalar*) ctx;
304: const PetscScalar x0 = vertices[0];
305: const PetscScalar y0 = vertices[1];
306: const PetscScalar x1 = vertices[2];
307: const PetscScalar y1 = vertices[3];
308: const PetscScalar x2 = vertices[4];
309: const PetscScalar y2 = vertices[5];
310: const PetscScalar x3 = vertices[6];
311: const PetscScalar y3 = vertices[7];
312: const PetscScalar f_1 = x1 - x0;
313: const PetscScalar g_1 = y1 - y0;
314: const PetscScalar f_3 = x3 - x0;
315: const PetscScalar g_3 = y3 - y0;
316: const PetscScalar f_01 = x2 - x1 - x3 + x0;
317: const PetscScalar g_01 = y2 - y1 - y3 + y0;
318: PetscScalar *ref, *real;
319: PetscErrorCode ierr;
322: VecGetArray(Xref, &ref);
323: VecGetArray(Xreal, &real);
324: {
325: const PetscScalar p0 = ref[0];
326: const PetscScalar p1 = ref[1];
328: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
329: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
330: }
331: PetscLogFlops(28);
332: VecRestoreArray(Xref, &ref);
333: VecRestoreArray(Xreal, &real);
334: return(0);
335: }
337: #include <petsc-private/dmimpl.h>
340: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
341: {
342: const PetscScalar *vertices = (const PetscScalar*) ctx;
343: const PetscScalar x0 = vertices[0];
344: const PetscScalar y0 = vertices[1];
345: const PetscScalar x1 = vertices[2];
346: const PetscScalar y1 = vertices[3];
347: const PetscScalar x2 = vertices[4];
348: const PetscScalar y2 = vertices[5];
349: const PetscScalar x3 = vertices[6];
350: const PetscScalar y3 = vertices[7];
351: const PetscScalar f_01 = x2 - x1 - x3 + x0;
352: const PetscScalar g_01 = y2 - y1 - y3 + y0;
353: PetscScalar *ref;
354: PetscErrorCode ierr;
357: VecGetArray(Xref, &ref);
358: {
359: const PetscScalar x = ref[0];
360: const PetscScalar y = ref[1];
361: const PetscInt rows[2] = {0, 1};
362: PetscScalar values[4];
364: values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
365: values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
366: MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
367: }
368: PetscLogFlops(30);
369: VecRestoreArray(Xref, &ref);
370: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
371: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
372: return(0);
373: }
377: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
378: {
379: DM dmCoord;
380: SNES snes;
381: KSP ksp;
382: PC pc;
383: Vec coordsLocal, r, ref, real;
384: Mat J;
385: PetscScalar *a, *coords;
386: PetscInt p;
390: DMGetCoordinatesLocal(dm, &coordsLocal);
391: DMGetCoordinateDM(dm, &dmCoord);
392: SNESCreate(PETSC_COMM_SELF, &snes);
393: SNESSetOptionsPrefix(snes, "quad_interp_");
394: VecCreate(PETSC_COMM_SELF, &r);
395: VecSetSizes(r, 2, 2);
396: VecSetType(r,dm->vectype);
397: VecDuplicate(r, &ref);
398: VecDuplicate(r, &real);
399: MatCreate(PETSC_COMM_SELF, &J);
400: MatSetSizes(J, 2, 2, 2, 2);
401: MatSetType(J, MATSEQDENSE);
402: MatSetUp(J);
403: SNESSetFunction(snes, r, QuadMap_Private, NULL);
404: SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
405: SNESGetKSP(snes, &ksp);
406: KSPGetPC(ksp, &pc);
407: PCSetType(pc, PCLU);
408: SNESSetFromOptions(snes);
410: VecGetArray(ctx->coords, &coords);
411: VecGetArray(v, &a);
412: for (p = 0; p < ctx->n; ++p) {
413: PetscScalar *x = NULL, *vertices = NULL;
414: PetscScalar *xi;
415: PetscReal xir[2];
416: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
418: /* Can make this do all points at once */
419: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
420: if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);
421: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
422: if (4*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 4*ctx->dof);
423: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
424: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
425: VecGetArray(real, &xi);
426: xi[0] = coords[p*ctx->dim+0];
427: xi[1] = coords[p*ctx->dim+1];
428: VecRestoreArray(real, &xi);
429: SNESSolve(snes, real, ref);
430: VecGetArray(ref, &xi);
431: xir[0] = PetscRealPart(xi[0]);
432: xir[1] = PetscRealPart(xi[1]);
433: 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];
435: VecRestoreArray(ref, &xi);
436: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
437: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
438: }
439: VecRestoreArray(v, &a);
440: VecRestoreArray(ctx->coords, &coords);
442: SNESDestroy(&snes);
443: VecDestroy(&r);
444: VecDestroy(&ref);
445: VecDestroy(&real);
446: MatDestroy(&J);
447: return(0);
448: }
452: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
453: {
454: const PetscScalar *vertices = (const PetscScalar*) ctx;
455: const PetscScalar x0 = vertices[0];
456: const PetscScalar y0 = vertices[1];
457: const PetscScalar z0 = vertices[2];
458: const PetscScalar x1 = vertices[9];
459: const PetscScalar y1 = vertices[10];
460: const PetscScalar z1 = vertices[11];
461: const PetscScalar x2 = vertices[6];
462: const PetscScalar y2 = vertices[7];
463: const PetscScalar z2 = vertices[8];
464: const PetscScalar x3 = vertices[3];
465: const PetscScalar y3 = vertices[4];
466: const PetscScalar z3 = vertices[5];
467: const PetscScalar x4 = vertices[12];
468: const PetscScalar y4 = vertices[13];
469: const PetscScalar z4 = vertices[14];
470: const PetscScalar x5 = vertices[15];
471: const PetscScalar y5 = vertices[16];
472: const PetscScalar z5 = vertices[17];
473: const PetscScalar x6 = vertices[18];
474: const PetscScalar y6 = vertices[19];
475: const PetscScalar z6 = vertices[20];
476: const PetscScalar x7 = vertices[21];
477: const PetscScalar y7 = vertices[22];
478: const PetscScalar z7 = vertices[23];
479: const PetscScalar f_1 = x1 - x0;
480: const PetscScalar g_1 = y1 - y0;
481: const PetscScalar h_1 = z1 - z0;
482: const PetscScalar f_3 = x3 - x0;
483: const PetscScalar g_3 = y3 - y0;
484: const PetscScalar h_3 = z3 - z0;
485: const PetscScalar f_4 = x4 - x0;
486: const PetscScalar g_4 = y4 - y0;
487: const PetscScalar h_4 = z4 - z0;
488: const PetscScalar f_01 = x2 - x1 - x3 + x0;
489: const PetscScalar g_01 = y2 - y1 - y3 + y0;
490: const PetscScalar h_01 = z2 - z1 - z3 + z0;
491: const PetscScalar f_12 = x7 - x3 - x4 + x0;
492: const PetscScalar g_12 = y7 - y3 - y4 + y0;
493: const PetscScalar h_12 = z7 - z3 - z4 + z0;
494: const PetscScalar f_02 = x5 - x1 - x4 + x0;
495: const PetscScalar g_02 = y5 - y1 - y4 + y0;
496: const PetscScalar h_02 = z5 - z1 - z4 + z0;
497: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
498: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
499: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
500: PetscScalar *ref, *real;
501: PetscErrorCode ierr;
504: VecGetArray(Xref, &ref);
505: VecGetArray(Xreal, &real);
506: {
507: const PetscScalar p0 = ref[0];
508: const PetscScalar p1 = ref[1];
509: const PetscScalar p2 = ref[2];
511: 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;
512: 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;
513: 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;
514: }
515: PetscLogFlops(114);
516: VecRestoreArray(Xref, &ref);
517: VecRestoreArray(Xreal, &real);
518: return(0);
519: }
523: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
524: {
525: const PetscScalar *vertices = (const PetscScalar*) ctx;
526: const PetscScalar x0 = vertices[0];
527: const PetscScalar y0 = vertices[1];
528: const PetscScalar z0 = vertices[2];
529: const PetscScalar x1 = vertices[9];
530: const PetscScalar y1 = vertices[10];
531: const PetscScalar z1 = vertices[11];
532: const PetscScalar x2 = vertices[6];
533: const PetscScalar y2 = vertices[7];
534: const PetscScalar z2 = vertices[8];
535: const PetscScalar x3 = vertices[3];
536: const PetscScalar y3 = vertices[4];
537: const PetscScalar z3 = vertices[5];
538: const PetscScalar x4 = vertices[12];
539: const PetscScalar y4 = vertices[13];
540: const PetscScalar z4 = vertices[14];
541: const PetscScalar x5 = vertices[15];
542: const PetscScalar y5 = vertices[16];
543: const PetscScalar z5 = vertices[17];
544: const PetscScalar x6 = vertices[18];
545: const PetscScalar y6 = vertices[19];
546: const PetscScalar z6 = vertices[20];
547: const PetscScalar x7 = vertices[21];
548: const PetscScalar y7 = vertices[22];
549: const PetscScalar z7 = vertices[23];
550: const PetscScalar f_xy = x2 - x1 - x3 + x0;
551: const PetscScalar g_xy = y2 - y1 - y3 + y0;
552: const PetscScalar h_xy = z2 - z1 - z3 + z0;
553: const PetscScalar f_yz = x7 - x3 - x4 + x0;
554: const PetscScalar g_yz = y7 - y3 - y4 + y0;
555: const PetscScalar h_yz = z7 - z3 - z4 + z0;
556: const PetscScalar f_xz = x5 - x1 - x4 + x0;
557: const PetscScalar g_xz = y5 - y1 - y4 + y0;
558: const PetscScalar h_xz = z5 - z1 - z4 + z0;
559: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
560: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
561: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
562: PetscScalar *ref;
563: PetscErrorCode ierr;
566: VecGetArray(Xref, &ref);
567: {
568: const PetscScalar x = ref[0];
569: const PetscScalar y = ref[1];
570: const PetscScalar z = ref[2];
571: const PetscInt rows[3] = {0, 1, 2};
572: PetscScalar values[9];
574: values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
575: values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
576: values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
577: values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
578: values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
579: values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
580: values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
581: values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
582: values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
584: MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
585: }
586: PetscLogFlops(152);
587: VecRestoreArray(Xref, &ref);
588: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
589: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
590: return(0);
591: }
595: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
596: {
597: DM dmCoord;
598: SNES snes;
599: KSP ksp;
600: PC pc;
601: Vec coordsLocal, r, ref, real;
602: Mat J;
603: PetscScalar *a, *coords;
604: PetscInt p;
608: DMGetCoordinatesLocal(dm, &coordsLocal);
609: DMGetCoordinateDM(dm, &dmCoord);
610: SNESCreate(PETSC_COMM_SELF, &snes);
611: SNESSetOptionsPrefix(snes, "hex_interp_");
612: VecCreate(PETSC_COMM_SELF, &r);
613: VecSetSizes(r, 3, 3);
614: VecSetType(r,dm->vectype);
615: VecDuplicate(r, &ref);
616: VecDuplicate(r, &real);
617: MatCreate(PETSC_COMM_SELF, &J);
618: MatSetSizes(J, 3, 3, 3, 3);
619: MatSetType(J, MATSEQDENSE);
620: MatSetUp(J);
621: SNESSetFunction(snes, r, HexMap_Private, NULL);
622: SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
623: SNESGetKSP(snes, &ksp);
624: KSPGetPC(ksp, &pc);
625: PCSetType(pc, PCLU);
626: SNESSetFromOptions(snes);
628: VecGetArray(ctx->coords, &coords);
629: VecGetArray(v, &a);
630: for (p = 0; p < ctx->n; ++p) {
631: PetscScalar *x = NULL, *vertices = NULL;
632: PetscScalar *xi;
633: PetscReal xir[3];
634: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
636: /* Can make this do all points at once */
637: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
638: if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);
639: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
640: if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);
641: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
642: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
643: VecGetArray(real, &xi);
644: xi[0] = coords[p*ctx->dim+0];
645: xi[1] = coords[p*ctx->dim+1];
646: xi[2] = coords[p*ctx->dim+2];
647: VecRestoreArray(real, &xi);
648: SNESSolve(snes, real, ref);
649: VecGetArray(ref, &xi);
650: xir[0] = PetscRealPart(xi[0]);
651: xir[1] = PetscRealPart(xi[1]);
652: xir[2] = PetscRealPart(xi[2]);
653: for (comp = 0; comp < ctx->dof; ++comp) {
654: a[p*ctx->dof+comp] =
655: x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
656: x[3*ctx->dof+comp]* xir[0]*(1-xir[1])*(1-xir[2]) +
657: x[2*ctx->dof+comp]* xir[0]* xir[1]*(1-xir[2]) +
658: x[1*ctx->dof+comp]*(1-xir[0])* xir[1]*(1-xir[2]) +
659: x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])* xir[2] +
660: x[5*ctx->dof+comp]* xir[0]*(1-xir[1])* xir[2] +
661: x[6*ctx->dof+comp]* xir[0]* xir[1]* xir[2] +
662: x[7*ctx->dof+comp]*(1-xir[0])* xir[1]* xir[2];
663: }
664: VecRestoreArray(ref, &xi);
665: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
666: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
667: }
668: VecRestoreArray(v, &a);
669: VecRestoreArray(ctx->coords, &coords);
671: SNESDestroy(&snes);
672: VecDestroy(&r);
673: VecDestroy(&ref);
674: VecDestroy(&real);
675: MatDestroy(&J);
676: return(0);
677: }
681: /*
682: Input Parameters:
683: + ctx - The DMInterpolationInfo context
684: . dm - The DM
685: - x - The local vector containing the field to be interpolated
687: Output Parameters:
688: . v - The vector containing the interpolated values
689: */
690: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
691: {
692: PetscInt dim, coneSize, n;
699: VecGetLocalSize(v, &n);
700: 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);
701: if (n) {
702: DMPlexGetDimension(dm, &dim);
703: DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);
704: if (dim == 2) {
705: if (coneSize == 3) {
706: DMInterpolate_Triangle_Private(ctx, dm, x, v);
707: } else if (coneSize == 4) {
708: DMInterpolate_Quad_Private(ctx, dm, x, v);
709: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
710: } else if (dim == 3) {
711: if (coneSize == 4) {
712: DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
713: } else {
714: DMInterpolate_Hex_Private(ctx, dm, x, v);
715: }
716: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
717: }
718: return(0);
719: }
723: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
724: {
729: VecDestroy(&(*ctx)->coords);
730: PetscFree((*ctx)->points);
731: PetscFree((*ctx)->cells);
732: PetscFree(*ctx);
733: *ctx = NULL;
734: return(0);
735: }
739: /*@C
740: SNESMonitorFields - Monitors the residual for each field separately
742: Collective on SNES
744: Input Parameters:
745: + snes - the SNES context
746: . its - iteration number
747: . fgnorm - 2-norm of residual
748: - dummy - unused context
750: Notes:
751: This routine prints the residual norm at each iteration.
753: Level: intermediate
755: .keywords: SNES, nonlinear, default, monitor, norm
756: .seealso: SNESMonitorSet(), SNESMonitorDefault()
757: @*/
758: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
759: {
760: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) snes));
761: Vec res;
762: DM dm;
763: PetscSection s;
764: const PetscScalar *r;
765: PetscReal *lnorms, *norms;
766: PetscInt numFields, f, pStart, pEnd, p;
767: PetscErrorCode ierr;
770: SNESGetFunction(snes, &res, 0, 0);
771: SNESGetDM(snes, &dm);
772: DMGetDefaultSection(dm, &s);
773: PetscSectionGetNumFields(s, &numFields);
774: PetscSectionGetChart(s, &pStart, &pEnd);
775: PetscCalloc2(numFields, &lnorms, numFields, &norms);
776: VecGetArrayRead(res, &r);
777: for (p = pStart; p < pEnd; ++p) {
778: for (f = 0; f < numFields; ++f) {
779: PetscInt fdof, foff, d;
781: PetscSectionGetFieldDof(s, p, f, &fdof);
782: PetscSectionGetFieldOffset(s, p, f, &foff);
783: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
784: }
785: }
786: VecRestoreArrayRead(res, &r);
787: MPI_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject) dm));
788: PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
789: PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
790: for (f = 0; f < numFields; ++f) {
791: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
792: PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
793: }
794: PetscViewerASCIIPrintf(viewer, "]\n");
795: PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
796: PetscFree2(lnorms, norms);
797: return(0);
798: }