Actual source code: dmmeshsnes.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/meshimpl.h> /*I "petscdmmesh.h" I*/
2: #include <petscsnes.h> /*I "petscsnes.h" I*/
6: PetscErrorCode DMMeshInterpolationCreate(DM dm, DMMeshInterpolationInfo *ctx)
7: {
13: PetscMalloc(sizeof(struct _DMMeshInterpolationInfo), ctx);
15: (*ctx)->dim = -1;
16: (*ctx)->nInput = 0;
17: (*ctx)->points = NULL;
18: (*ctx)->cells = NULL;
19: (*ctx)->n = -1;
20: (*ctx)->coords = NULL;
21: return(0);
22: }
26: PetscErrorCode DMMeshInterpolationSetDim(DM dm, PetscInt dim, DMMeshInterpolationInfo ctx)
27: {
30: if ((dim < 1) || (dim > 3)) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);
31: ctx->dim = dim;
32: return(0);
33: }
37: PetscErrorCode DMMeshInterpolationGetDim(DM dm, PetscInt *dim, DMMeshInterpolationInfo ctx)
38: {
42: *dim = ctx->dim;
43: return(0);
44: }
48: PetscErrorCode DMMeshInterpolationSetDof(DM dm, PetscInt dof, DMMeshInterpolationInfo ctx)
49: {
52: if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);
53: ctx->dof = dof;
54: return(0);
55: }
59: PetscErrorCode DMMeshInterpolationGetDof(DM dm, PetscInt *dof, DMMeshInterpolationInfo ctx)
60: {
64: *dof = ctx->dof;
65: return(0);
66: }
70: PetscErrorCode DMMeshInterpolationAddPoints(DM dm, PetscInt n, PetscReal points[], DMMeshInterpolationInfo ctx)
71: {
76: if (ctx->dim < 0) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
77: if (ctx->points) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
78: ctx->nInput = n;
79: PetscMalloc(n*ctx->dim * sizeof(PetscReal), &ctx->points);
80: PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));
81: return(0);
82: }
86: PetscErrorCode DMMeshInterpolationSetUp(DM dm, DMMeshInterpolationInfo ctx, PetscBool redundantPoints)
87: {
88: ALE::Obj<PETSC_MESH_TYPE> m;
89: MPI_Comm comm;
90: PetscScalar *a;
91: PetscInt p, q, i;
92: PetscMPIInt rank, size;
93: PetscErrorCode ierr;
97: PetscObjectGetComm((PetscObject)dm,&comm);
98: MPI_Comm_size(comm, &size);
99: MPI_Comm_rank(comm, &rank);
100: DMMeshGetMesh(dm, m);
101: if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
103: /* Locate points */
104: PetscLayout layout;
105: PetscReal *globalPoints;
106: const PetscInt *ranges;
107: PetscMPIInt *counts, *displs;
108: PetscInt *foundCells;
109: PetscMPIInt *foundProcs, *globalProcs;
110: PetscInt n = ctx->nInput, N;
112: if (!redundantPoints) {
113: PetscLayoutCreate(comm, &layout);
114: PetscLayoutSetBlockSize(layout, 1);
115: PetscLayoutSetLocalSize(layout, n);
116: PetscLayoutSetUp(layout);
117: PetscLayoutGetSize(layout, &N);
118: /* Communicate all points to all processes */
119: PetscMalloc3(N*ctx->dim,PetscReal,&globalPoints,size,PetscMPIInt,&counts,size,PetscMPIInt,&displs);
120: PetscLayoutGetRanges(layout, &ranges);
121: for (p = 0; p < size; ++p) {
122: counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
123: displs[p] = ranges[p]*ctx->dim;
124: }
125: MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
126: } else {
127: N = n;
129: globalPoints = ctx->points;
130: }
131: PetscMalloc3(N,PetscInt,&foundCells,N,PetscMPIInt,&foundProcs,N,PetscMPIInt,&globalProcs);
132: for (p = 0; p < N; ++p) {
133: foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]);
134: if (foundCells[p] >= 0) foundProcs[p] = rank;
135: else foundProcs[p] = size;
136: }
137: /* Let the lowest rank process own each point */
138: MPI_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
139: ctx->n = 0;
140: for (p = 0; p < N; ++p) {
141: 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);
142: else if (globalProcs[p] == rank) ctx->n++;
143: }
144: /* Create coordinates vector and array of owned cells */
145: PetscMalloc(ctx->n * sizeof(PetscInt), &ctx->cells);
146: VecCreate(comm, &ctx->coords);
147: VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
148: VecSetBlockSize(ctx->coords, ctx->dim);
149: VecSetFromOptions(ctx->coords);
150: VecGetArray(ctx->coords, &a);
151: for (p = 0, q = 0, i = 0; p < N; ++p) {
152: if (globalProcs[p] == rank) {
153: PetscInt d;
155: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
156: ctx->cells[q++] = foundCells[p];
157: }
158: }
159: VecRestoreArray(ctx->coords, &a);
160: PetscFree3(foundCells,foundProcs,globalProcs);
161: if (!redundantPoints) {
162: PetscFree3(globalPoints,counts,displs);
163: PetscLayoutDestroy(&layout);
164: }
165: return(0);
166: }
170: PetscErrorCode DMMeshInterpolationGetCoordinates(DM dm, Vec *coordinates, DMMeshInterpolationInfo ctx)
171: {
175: if (!ctx->coords) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
176: *coordinates = ctx->coords;
177: return(0);
178: }
182: PetscErrorCode DMMeshInterpolationGetVector(DM dm, Vec *v, DMMeshInterpolationInfo ctx)
183: {
189: if (!ctx->coords) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
190: VecCreate(((PetscObject) dm)->comm, v);
191: VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
192: VecSetBlockSize(*v, ctx->dof);
193: VecSetFromOptions(*v);
194: return(0);
195: }
199: PetscErrorCode DMMeshInterpolationRestoreVector(DM dm, Vec *v, DMMeshInterpolationInfo ctx)
200: {
206: if (!ctx->coords) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
207: VecDestroy(v);
208: return(0);
209: }
213: PETSC_STATIC_INLINE PetscErrorCode DMMeshInterpolate_Simplex_Private(DM dm, SectionReal x, Vec v, DMMeshInterpolationInfo ctx)
214: {
215: #if defined(PETSC_HAVE_SIEVE)
216: ALE::Obj<PETSC_MESH_TYPE> m;
217: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
218: PetscInt p;
219: PetscErrorCode ierr;
222: DMMeshGetMesh(dm, m);
223: SectionRealGetSection(x, s);
224: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = m->getRealSection("coordinates");
225: PetscReal *v0, *J, *invJ, detJ;
226: PetscScalar *a, *coords;
228: PetscMalloc3(ctx->dim,PetscReal,&v0,ctx->dim*ctx->dim,PetscReal,&J,ctx->dim*ctx->dim,PetscReal,&invJ);
229: VecGetArray(ctx->coords, &coords);
230: VecGetArray(v, &a);
231: for (p = 0; p < ctx->n; ++p) {
232: PetscInt e = ctx->cells[p];
233: PetscReal xi[4];
234: PetscInt d, f, comp;
236: if ((ctx->dim+1)*ctx->dof != m->sizeWithBC(s, e)) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_ARG_SIZ, "Invalid restrict size %d should be %d", m->sizeWithBC(s, e), (ctx->dim+1)*ctx->dof);
237: m->computeElementGeometry(coordinates, e, v0, J, invJ, detJ);
238: const PetscScalar *c = m->restrictClosure(s, e); /* Must come after geom, since it uses closure temp space*/
239: for (comp = 0; comp < ctx->dof; ++comp) {
240: a[p*ctx->dof+comp] = c[0*ctx->dof+comp];
241: }
242: for (d = 0; d < ctx->dim; ++d) {
243: xi[d] = 0.0;
244: for (f = 0; f < ctx->dim; ++f) {
245: xi[d] += invJ[d*ctx->dim+f]*0.5*(coords[p*ctx->dim+f] - v0[f]);
246: }
247: for (comp = 0; comp < ctx->dof; ++comp) {
248: a[p*ctx->dof+comp] += (c[(d+1)*ctx->dof+comp] - c[0*ctx->dof+comp])*xi[d];
249: }
250: }
251: }
252: VecRestoreArray(v, &a);
253: VecRestoreArray(ctx->coords, &coords);
254: PetscFree3(v0, J, invJ);
255: return(0);
256: #else
257: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Interpolation only work with DMMesh currently.");
258: #endif
259: }
263: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
264: {
265: const PetscScalar *vertices = (const PetscScalar*) ctx;
266: const PetscScalar x0 = vertices[0];
267: const PetscScalar y0 = vertices[1];
268: const PetscScalar x1 = vertices[2];
269: const PetscScalar y1 = vertices[3];
270: const PetscScalar x2 = vertices[4];
271: const PetscScalar y2 = vertices[5];
272: const PetscScalar x3 = vertices[6];
273: const PetscScalar y3 = vertices[7];
274: const PetscScalar f_1 = x1 - x0;
275: const PetscScalar g_1 = y1 - y0;
276: const PetscScalar f_3 = x3 - x0;
277: const PetscScalar g_3 = y3 - y0;
278: const PetscScalar f_01 = x2 - x1 - x3 + x0;
279: const PetscScalar g_01 = y2 - y1 - y3 + y0;
280: PetscScalar *ref, *real;
281: PetscErrorCode ierr;
284: VecGetArray(Xref, &ref);
285: VecGetArray(Xreal, &real);
286: {
287: const PetscScalar p0 = ref[0];
288: const PetscScalar p1 = ref[1];
290: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
291: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
292: }
293: PetscLogFlops(28);
294: VecRestoreArray(Xref, &ref);
295: VecRestoreArray(Xreal, &real);
296: return(0);
297: }
301: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat *J, Mat *M, MatStructure *flag, 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_01 = x2 - x1 - x3 + x0;
313: const PetscScalar g_01 = y2 - y1 - y3 + y0;
314: PetscScalar *ref;
315: PetscErrorCode ierr;
318: VecGetArray(Xref, &ref);
319: {
320: const PetscScalar x = ref[0];
321: const PetscScalar y = ref[1];
322: const PetscInt rows[2] = {0, 1};
323: const PetscScalar values[4] = {(x1 - x0 + f_01*y) * 0.5, (x3 - x0 + f_01*x) * 0.5,
324: (y1 - y0 + g_01*y) * 0.5, (y3 - y0 + g_01*x) * 0.5};
325: MatSetValues(*J, 2, rows, 2, rows, values, INSERT_VALUES);
326: }
327: PetscLogFlops(30);
328: VecRestoreArray(Xref, &ref);
329: MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);
330: MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY);
331: return(0);
332: }
336: PETSC_STATIC_INLINE PetscErrorCode DMMeshInterpolate_Quad_Private(DM dm, SectionReal x, Vec v, DMMeshInterpolationInfo ctx)
337: {
338: #if defined(PETSC_HAVE_SIEVE)
339: SNES snes;
340: KSP ksp;
341: PC pc;
342: Vec r, ref, real;
343: Mat J;
344: PetscScalar vertices[8];
346: ALE::Obj<PETSC_MESH_TYPE> m;
347: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
348: PetscInt p;
349: PetscErrorCode ierr;
352: SNESCreate(PETSC_COMM_SELF, &snes);
353: SNESSetOptionsPrefix(snes, "quad_interp_");
354: VecCreate(PETSC_COMM_SELF, &r);
355: VecSetSizes(r, 2, 2);
356: VecSetFromOptions(r);
357: VecDuplicate(r, &ref);
358: VecDuplicate(r, &real);
359: MatCreate(PETSC_COMM_SELF, &J);
360: MatSetSizes(J, 2, 2, 2, 2);
361: MatSetType(J, MATSEQDENSE);
362: MatSetUp(J);
363: SNESSetFunction(snes, r, QuadMap_Private, vertices);
364: SNESSetJacobian(snes, J, J, QuadJacobian_Private, vertices);
365: SNESGetKSP(snes, &ksp);
366: KSPGetPC(ksp, &pc);
367: PCSetType(pc, PCLU);
368: SNESSetFromOptions(snes);
370: DMMeshGetMesh(dm, m);
371: SectionRealGetSection(x, s);
372: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = m->getRealSection("coordinates");
373: PetscScalar *a, *coords;
375: VecGetArray(ctx->coords, &coords);
376: VecGetArray(v, &a);
377: for (p = 0; p < ctx->n; ++p) {
378: PetscScalar *xi;
379: PetscInt e = ctx->cells[p], comp;
381: if (4*ctx->dof != m->sizeWithBC(s, e)) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_ARG_SIZ, "Invalid restrict size %d should be %d", m->sizeWithBC(s, e), 4*ctx->dof);
382: /* Can make this do all points at once */
383: {
384: const PetscReal *v = m->restrictClosure(coordinates, e);
385: for (PetscInt i = 0; i < 8; ++i) vertices[i] = v[i];
386: }
387: const PetscScalar *c = m->restrictClosure(s, e); /* Must come after geom, since it uses closure temp space*/
388: VecGetArray(real, &xi);
389: xi[0] = coords[p*ctx->dim+0];
390: xi[1] = coords[p*ctx->dim+1];
391: VecRestoreArray(real, &xi);
392: SNESSolve(snes, real, ref);
393: VecGetArray(ref, &xi);
394: for (comp = 0; comp < ctx->dof; ++comp) {
395: a[p*ctx->dof+comp] = c[0*ctx->dof+comp]*(1 - xi[0])*(1 - xi[1]) + c[1*ctx->dof+comp]*xi[0]*(1 - xi[1]) + c[2*ctx->dof+comp]*xi[0]*xi[1] + c[3*ctx->dof+comp]*(1 - xi[0])*xi[1];
396: }
397: VecRestoreArray(ref, &xi);
398: }
399: VecRestoreArray(v, &a);
400: VecRestoreArray(ctx->coords, &coords);
402: SNESDestroy(&snes);
403: VecDestroy(&r);
404: VecDestroy(&ref);
405: VecDestroy(&real);
406: MatDestroy(&J);
407: return(0);
408: #else
409: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Interpolation only work with DMMesh currently.");
410: #endif
411: }
415: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
416: {
417: const PetscScalar *vertices = (const PetscScalar*) ctx;
418: const PetscScalar x0 = vertices[0];
419: const PetscScalar y0 = vertices[1];
420: const PetscScalar z0 = vertices[2];
421: const PetscScalar x1 = vertices[3];
422: const PetscScalar y1 = vertices[4];
423: const PetscScalar z1 = vertices[5];
424: const PetscScalar x2 = vertices[6];
425: const PetscScalar y2 = vertices[7];
426: const PetscScalar z2 = vertices[8];
427: const PetscScalar x3 = vertices[9];
428: const PetscScalar y3 = vertices[10];
429: const PetscScalar z3 = vertices[11];
430: const PetscScalar x4 = vertices[12];
431: const PetscScalar y4 = vertices[13];
432: const PetscScalar z4 = vertices[14];
433: const PetscScalar x5 = vertices[15];
434: const PetscScalar y5 = vertices[16];
435: const PetscScalar z5 = vertices[17];
436: const PetscScalar x6 = vertices[18];
437: const PetscScalar y6 = vertices[19];
438: const PetscScalar z6 = vertices[20];
439: const PetscScalar x7 = vertices[21];
440: const PetscScalar y7 = vertices[22];
441: const PetscScalar z7 = vertices[23];
442: const PetscScalar f_1 = x1 - x0;
443: const PetscScalar g_1 = y1 - y0;
444: const PetscScalar h_1 = z1 - z0;
445: const PetscScalar f_3 = x3 - x0;
446: const PetscScalar g_3 = y3 - y0;
447: const PetscScalar h_3 = z3 - z0;
448: const PetscScalar f_4 = x4 - x0;
449: const PetscScalar g_4 = y4 - y0;
450: const PetscScalar h_4 = z4 - z0;
451: const PetscScalar f_01 = x2 - x1 - x3 + x0;
452: const PetscScalar g_01 = y2 - y1 - y3 + y0;
453: const PetscScalar h_01 = z2 - z1 - z3 + z0;
454: const PetscScalar f_12 = x7 - x3 - x4 + x0;
455: const PetscScalar g_12 = y7 - y3 - y4 + y0;
456: const PetscScalar h_12 = z7 - z3 - z4 + z0;
457: const PetscScalar f_02 = x5 - x1 - x4 + x0;
458: const PetscScalar g_02 = y5 - y1 - y4 + y0;
459: const PetscScalar h_02 = z5 - z1 - z4 + z0;
460: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
461: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
462: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
463: PetscScalar *ref, *real;
464: PetscErrorCode ierr;
467: VecGetArray(Xref, &ref);
468: VecGetArray(Xreal, &real);
469: {
470: const PetscScalar p0 = ref[0];
471: const PetscScalar p1 = ref[1];
472: const PetscScalar p2 = ref[2];
474: 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;
475: 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;
476: 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;
477: }
478: PetscLogFlops(114);
479: VecRestoreArray(Xref, &ref);
480: VecRestoreArray(Xreal, &real);
481: return(0);
482: }
486: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat *J, Mat *M, MatStructure *flag, void *ctx)
487: {
488: const PetscScalar *vertices = (const PetscScalar*) ctx;
489: const PetscScalar x0 = vertices[0];
490: const PetscScalar y0 = vertices[1];
491: const PetscScalar z0 = vertices[2];
492: const PetscScalar x1 = vertices[3];
493: const PetscScalar y1 = vertices[4];
494: const PetscScalar z1 = vertices[5];
495: const PetscScalar x2 = vertices[6];
496: const PetscScalar y2 = vertices[7];
497: const PetscScalar z2 = vertices[8];
498: const PetscScalar x3 = vertices[9];
499: const PetscScalar y3 = vertices[10];
500: const PetscScalar z3 = vertices[11];
501: const PetscScalar x4 = vertices[12];
502: const PetscScalar y4 = vertices[13];
503: const PetscScalar z4 = vertices[14];
504: const PetscScalar x5 = vertices[15];
505: const PetscScalar y5 = vertices[16];
506: const PetscScalar z5 = vertices[17];
507: const PetscScalar x6 = vertices[18];
508: const PetscScalar y6 = vertices[19];
509: const PetscScalar z6 = vertices[20];
510: const PetscScalar x7 = vertices[21];
511: const PetscScalar y7 = vertices[22];
512: const PetscScalar z7 = vertices[23];
513: const PetscScalar f_xy = x2 - x1 - x3 + x0;
514: const PetscScalar g_xy = y2 - y1 - y3 + y0;
515: const PetscScalar h_xy = z2 - z1 - z3 + z0;
516: const PetscScalar f_yz = x7 - x3 - x4 + x0;
517: const PetscScalar g_yz = y7 - y3 - y4 + y0;
518: const PetscScalar h_yz = z7 - z3 - z4 + z0;
519: const PetscScalar f_xz = x5 - x1 - x4 + x0;
520: const PetscScalar g_xz = y5 - y1 - y4 + y0;
521: const PetscScalar h_xz = z5 - z1 - z4 + z0;
522: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
523: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
524: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
525: PetscScalar *ref;
526: PetscErrorCode ierr;
529: VecGetArray(Xref, &ref);
530: {
531: const PetscScalar x = ref[0];
532: const PetscScalar y = ref[1];
533: const PetscScalar z = ref[2];
534: const PetscInt rows[3] = {0, 1, 2};
535: const PetscScalar values[9] = {
536: (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0,
537: (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0,
538: (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0,
539: (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0,
540: (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0,
541: (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0,
542: (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0,
543: (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0,
544: (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0
545: };
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 DMMeshInterpolate_Hex_Private(DM dm, SectionReal x, Vec v, DMMeshInterpolationInfo ctx)
558: {
559: #if defined(PETSC_HAVE_SIEVE)
560: SNES snes;
561: KSP ksp;
562: PC pc;
563: Vec r, ref, real;
564: Mat J;
565: PetscScalar vertices[24];
567: ALE::Obj<PETSC_MESH_TYPE> m;
568: ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
569: PetscInt p;
570: PetscErrorCode ierr;
573: SNESCreate(PETSC_COMM_SELF, &snes);
574: SNESSetOptionsPrefix(snes, "hex_interp_");
575: VecCreate(PETSC_COMM_SELF, &r);
576: VecSetSizes(r, 3, 3);
577: VecSetFromOptions(r);
578: VecDuplicate(r, &ref);
579: VecDuplicate(r, &real);
580: MatCreate(PETSC_COMM_SELF, &J);
581: MatSetSizes(J, 3, 3, 3, 3);
582: MatSetType(J, MATSEQDENSE);
583: MatSetUp(J);
584: SNESSetFunction(snes, r, HexMap_Private, vertices);
585: SNESSetJacobian(snes, J, J, HexJacobian_Private, vertices);
586: SNESGetKSP(snes, &ksp);
587: KSPGetPC(ksp, &pc);
588: PCSetType(pc, PCLU);
589: SNESSetFromOptions(snes);
591: DMMeshGetMesh(dm, m);
592: SectionRealGetSection(x, s);
593: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = m->getRealSection("coordinates");
594: PetscScalar *a, *coords;
596: VecGetArray(ctx->coords, &coords);
597: VecGetArray(v, &a);
598: for (p = 0; p < ctx->n; ++p) {
599: PetscScalar *xi;
600: PetscInt e = ctx->cells[p], comp;
602: if (8*ctx->dof != m->sizeWithBC(s, e)) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_ARG_SIZ, "Invalid restrict size %d should be %d", m->sizeWithBC(s, e), 8*ctx->dof);
603: /* Can make this do all points at once */
604: {
605: const PetscReal *v = m->restrictClosure(coordinates, e);
606: for (PetscInt i = 0; i < 24; ++i) vertices[i] = v[i];
607: }
608: const PetscScalar *c = m->restrictClosure(s, e); /* Must come after geom, since it uses closure temp space*/
609: VecGetArray(real, &xi);
610: xi[0] = coords[p*ctx->dim+0];
611: xi[1] = coords[p*ctx->dim+1];
612: xi[2] = coords[p*ctx->dim+2];
613: VecRestoreArray(real, &xi);
614: SNESSolve(snes, real, ref);
615: VecGetArray(ref, &xi);
616: for (comp = 0; comp < ctx->dof; ++comp) {
617: a[p*ctx->dof+comp] =
618: c[0*ctx->dof+comp]*(1-xi[0])*(1-xi[1])*(1-xi[2]) +
619: c[1*ctx->dof+comp]* xi[0]*(1-xi[1])*(1-xi[2]) +
620: c[2*ctx->dof+comp]* xi[0]* xi[1]*(1-xi[2]) +
621: c[3*ctx->dof+comp]*(1-xi[0])* xi[1]*(1-xi[2]) +
622: c[4*ctx->dof+comp]*(1-xi[0])*(1-xi[1])* xi[2] +
623: c[5*ctx->dof+comp]* xi[0]*(1-xi[1])* xi[2] +
624: c[6*ctx->dof+comp]* xi[0]* xi[1]* xi[2] +
625: c[7*ctx->dof+comp]*(1-xi[0])* xi[1]* xi[2];
626: }
627: VecRestoreArray(ref, &xi);
628: }
629: VecRestoreArray(v, &a);
630: VecRestoreArray(ctx->coords, &coords);
632: SNESDestroy(&snes);
633: VecDestroy(&r);
634: VecDestroy(&ref);
635: VecDestroy(&real);
636: MatDestroy(&J);
637: return(0);
638: #else
639: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Interpolation only work with DMMesh currently.");
640: #endif
641: }
645: PetscErrorCode DMMeshInterpolationEvaluate(DM dm, SectionReal x, Vec v, DMMeshInterpolationInfo ctx)
646: {
647: PetscInt dim, coneSize, n;
654: VecGetLocalSize(v, &n);
655: if (n != ctx->n*ctx->dof) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %d should be %d", n, ctx->n*ctx->dof);
656: if (n) {
657: DMMeshGetDimension(dm, &dim);
658: DMMeshGetConeSize(dm, ctx->cells[0], &coneSize);
659: if (dim == 2) {
660: if (coneSize == 3) {
661: DMMeshInterpolate_Simplex_Private(dm, x, v, ctx);
662: } else if (coneSize == 4) {
663: DMMeshInterpolate_Quad_Private(dm, x, v, ctx);
664: } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
665: } else if (dim == 3) {
666: if (coneSize == 4) {
667: DMMeshInterpolate_Simplex_Private(dm, x, v, ctx);
668: } else {
669: DMMeshInterpolate_Hex_Private(dm, x, v, ctx);
670: }
671: } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
672: }
673: return(0);
674: }
678: PetscErrorCode DMMeshInterpolationDestroy(DM dm, DMMeshInterpolationInfo *ctx)
679: {
685: VecDestroy(&(*ctx)->coords);
686: PetscFree((*ctx)->points);
687: PetscFree((*ctx)->cells);
688: PetscFree(*ctx);
689: *ctx = NULL;
690: return(0);
691: }