Actual source code: plexfem.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petscsf.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petsc/private/petscfvimpl.h>
9: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
10: {
11: DM_Plex *mesh = (DM_Plex*) dm->data;
16: *scale = mesh->scale[unit];
17: return(0);
18: }
22: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
23: {
24: DM_Plex *mesh = (DM_Plex*) dm->data;
28: mesh->scale[unit] = scale;
29: return(0);
30: }
32: PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
33: {
34: switch (i) {
35: case 0:
36: switch (j) {
37: case 0: return 0;
38: case 1:
39: switch (k) {
40: case 0: return 0;
41: case 1: return 0;
42: case 2: return 1;
43: }
44: case 2:
45: switch (k) {
46: case 0: return 0;
47: case 1: return -1;
48: case 2: return 0;
49: }
50: }
51: case 1:
52: switch (j) {
53: case 0:
54: switch (k) {
55: case 0: return 0;
56: case 1: return 0;
57: case 2: return -1;
58: }
59: case 1: return 0;
60: case 2:
61: switch (k) {
62: case 0: return 1;
63: case 1: return 0;
64: case 2: return 0;
65: }
66: }
67: case 2:
68: switch (j) {
69: case 0:
70: switch (k) {
71: case 0: return 0;
72: case 1: return 1;
73: case 2: return 0;
74: }
75: case 1:
76: switch (k) {
77: case 0: return -1;
78: case 1: return 0;
79: case 2: return 0;
80: }
81: case 2: return 0;
82: }
83: }
84: return 0;
85: }
89: static PetscErrorCode DMPlexProjectRigidBody(PetscInt dim, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
90: {
91: PetscInt *ctxInt = (PetscInt *) ctx;
92: PetscInt dim2 = ctxInt[0];
93: PetscInt d = ctxInt[1];
94: PetscInt i, j, k = dim > 2 ? d - dim : d;
97: if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
98: for (i = 0; i < dim; i++) mode[i] = 0.;
99: if (d < dim) {
100: mode[d] = 1.;
101: } else {
102: for (i = 0; i < dim; i++) {
103: for (j = 0; j < dim; j++) {
104: mode[j] += epsilon(i, j, k)*X[i];
105: }
106: }
107: }
108: return(0);
109: }
113: /*@C
114: DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates
116: Collective on DM
118: Input Arguments:
119: . dm - the DM
121: Output Argument:
122: . sp - the null space
124: Note: This is necessary to take account of Dirichlet conditions on the displacements
126: Level: advanced
128: .seealso: MatNullSpaceCreate()
129: @*/
130: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
131: {
132: MPI_Comm comm;
133: Vec mode[6];
134: PetscSection section, globalSection;
135: PetscInt dim, dimEmbed, n, m, d, i, j;
139: PetscObjectGetComm((PetscObject)dm,&comm);
140: DMGetDimension(dm, &dim);
141: DMGetCoordinateDim(dm, &dimEmbed);
142: if (dim == 1) {
143: MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
144: return(0);
145: }
146: DMGetDefaultSection(dm, §ion);
147: DMGetDefaultGlobalSection(dm, &globalSection);
148: PetscSectionGetConstrainedStorageSize(globalSection, &n);
149: m = (dim*(dim+1))/2;
150: VecCreate(comm, &mode[0]);
151: VecSetSizes(mode[0], n, PETSC_DETERMINE);
152: VecSetUp(mode[0]);
153: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
154: for (d = 0; d < m; d++) {
155: PetscInt ctx[2];
156: void *voidctx = (void *) (&ctx[0]);
157: PetscErrorCode (*func)(PetscInt, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody;
159: ctx[0] = dimEmbed;
160: ctx[1] = d;
161: DMPlexProjectFunction(dm, &func, &voidctx, INSERT_VALUES, mode[d]);
162: }
163: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
164: /* Orthonormalize system */
165: for (i = dim; i < m; ++i) {
166: PetscScalar dots[6];
168: VecMDot(mode[i], i, mode, dots);
169: for (j = 0; j < i; ++j) dots[j] *= -1.0;
170: VecMAXPY(mode[i], i, dots, mode);
171: VecNormalize(mode[i], NULL);
172: }
173: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
174: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
175: return(0);
176: }
180: /*@
181: DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
182: are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
183: evaluating the dual space basis of that point. A basis function is associated with the point in its
184: transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
185: projection height, which is set with this function. By default, the maximum projection height is zero, which means
186: that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior
187: basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
189: Input Parameters:
190: + dm - the DMPlex object
191: - height - the maximum projection height >= 0
193: Level: advanced
195: .seealso: DMPlexGetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
196: @*/
197: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
198: {
199: DM_Plex *plex = (DM_Plex *) dm->data;
203: plex->maxProjectionHeight = height;
204: return(0);
205: }
209: /*@
210: DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
211: DMPlexProjectXXXLocal() functions.
213: Input Parameters:
214: . dm - the DMPlex object
216: Output Parameters:
217: . height - the maximum projection height
219: Level: intermediate
221: .seealso: DMPlexSetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
222: @*/
223: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
224: {
225: DM_Plex *plex = (DM_Plex *) dm->data;
229: *height = plex->maxProjectionHeight;
230: return(0);
231: }
235: PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
236: {
237: PetscDualSpace *sp, *cellsp;
238: PetscInt *numComp;
239: PetscSection section;
240: PetscScalar *values;
241: PetscBool *fieldActive;
242: PetscInt numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h;
243: PetscErrorCode ierr;
246: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
247: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
248: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
249: if (cEnd <= cStart) return(0);
250: DMGetDimension(dm, &dim);
251: DMGetCoordinateDim(dm, &dimEmbed);
252: DMGetDefaultSection(dm, §ion);
253: PetscSectionGetNumFields(section, &numFields);
254: PetscMalloc2(numFields,&sp,numFields,&numComp);
255: DMPlexGetMaxProjectionHeight(dm,&maxHeight);
256: if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
257: if (maxHeight > 0) {PetscMalloc1(numFields,&cellsp);}
258: else {cellsp = sp;}
259: for (h = 0; h <= maxHeight; h++) {
260: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
261: if (!h) {pStart = cStart; pEnd = cEnd;}
262: if (pEnd <= pStart) continue;
263: totDim = 0;
264: for (f = 0; f < numFields; ++f) {
265: PetscObject obj;
266: PetscClassId id;
268: DMGetField(dm, f, &obj);
269: PetscObjectGetClassId(obj, &id);
270: if (id == PETSCFE_CLASSID) {
271: PetscFE fe = (PetscFE) obj;
273: PetscFEGetNumComponents(fe, &numComp[f]);
274: if (!h) {
275: PetscFEGetDualSpace(fe, &cellsp[f]);
276: sp[f] = cellsp[f];
277: PetscObjectReference((PetscObject) sp[f]);
278: } else {
279: PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
280: if (!sp[f]) continue;
281: }
282: } else if (id == PETSCFV_CLASSID) {
283: PetscFV fv = (PetscFV) obj;
285: PetscFVGetNumComponents(fv, &numComp[f]);
286: PetscFVGetDualSpace(fv, &sp[f]);
287: PetscObjectReference((PetscObject) sp[f]);
288: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
289: PetscDualSpaceGetDimension(sp[f], &spDim);
290: totDim += spDim*numComp[f];
291: }
292: DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
293: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
294: if (!totDim) continue;
295: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
296: DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
297: for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
298: for (i = 0; i < numIds; ++i) {
299: IS pointIS;
300: const PetscInt *points;
301: PetscInt n, p;
303: DMLabelGetStratumIS(label, ids[i], &pointIS);
304: ISGetLocalSize(pointIS, &n);
305: ISGetIndices(pointIS, &points);
306: for (p = 0; p < n; ++p) {
307: const PetscInt point = points[p];
308: PetscFECellGeom geom;
310: if ((point < pStart) || (point >= pEnd)) continue;
311: DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);
312: geom.dim = dim - h;
313: geom.dimEmbed = dimEmbed;
314: for (f = 0, v = 0; f < numFields; ++f) {
315: void * const ctx = ctxs ? ctxs[f] : NULL;
317: if (!sp[f]) continue;
318: PetscDualSpaceGetDimension(sp[f], &spDim);
319: for (d = 0; d < spDim; ++d) {
320: if (funcs[f]) {
321: PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);
322: } else {
323: for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
324: }
325: v += numComp[f];
326: }
327: }
328: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);
329: }
330: ISRestoreIndices(pointIS, &points);
331: ISDestroy(&pointIS);
332: }
333: if (h) {
334: for (f = 0; f < numFields; ++f) {PetscDualSpaceDestroy(&sp[f]);}
335: }
336: }
337: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
338: DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
339: for (f = 0; f < numFields; ++f) {PetscDualSpaceDestroy(&sp[f]);}
340: PetscFree2(sp, numComp);
341: if (maxHeight > 0) {
342: PetscFree(cellsp);
343: }
344: return(0);
345: }
349: PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
350: {
351: PetscDualSpace *sp, *cellsp;
352: PetscInt *numComp;
353: PetscSection section;
354: PetscScalar *values;
355: PetscInt Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight;
356: PetscErrorCode ierr;
359: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
360: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
361: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
362: if (cEnd <= cStart) return(0);
363: DMGetDefaultSection(dm, §ion);
364: PetscSectionGetNumFields(section, &Nf);
365: PetscMalloc2(Nf, &sp, Nf, &numComp);
366: DMGetDimension(dm, &dim);
367: DMGetCoordinateDim(dm, &dimEmbed);
368: DMPlexGetMaxProjectionHeight(dm,&maxHeight);
369: if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
370: if (maxHeight > 0) {
371: PetscMalloc1(Nf,&cellsp);
372: }
373: else {
374: cellsp = sp;
375: }
376: for (h = 0; h <= maxHeight; h++) {
377: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
378: if (!h) {pStart = cStart; pEnd = cEnd;}
379: if (pEnd <= pStart) continue;
380: totDim = 0;
381: for (f = 0; f < Nf; ++f) {
382: PetscObject obj;
383: PetscClassId id;
385: DMGetField(dm, f, &obj);
386: PetscObjectGetClassId(obj, &id);
387: if (id == PETSCFE_CLASSID) {
388: PetscFE fe = (PetscFE) obj;
390: PetscFEGetNumComponents(fe, &numComp[f]);
391: if (!h) {
392: PetscFEGetDualSpace(fe, &cellsp[f]);
393: sp[f] = cellsp[f];
394: PetscObjectReference((PetscObject) sp[f]);
395: }
396: else {
397: PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
398: if (!sp[f]) {
399: continue;
400: }
401: }
402: } else if (id == PETSCFV_CLASSID) {
403: PetscFV fv = (PetscFV) obj;
405: if (h) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, "Projection height > 0 not supported for finite volume");
406: PetscFVGetNumComponents(fv, &numComp[f]);
407: PetscFVGetDualSpace(fv, &sp[f]);
408: PetscObjectReference((PetscObject) sp[f]);
409: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
410: PetscDualSpaceGetDimension(sp[f], &spDim);
411: totDim += spDim*numComp[f];
412: }
413: DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
414: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
415: if (!totDim) continue;
416: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
417: for (p = pStart; p < pEnd; ++p) {
418: PetscFECellGeom geom;
420: DMPlexComputeCellGeometryFEM(dm, p, NULL, geom.v0, geom.J, NULL, &geom.detJ);
421: geom.dim = dim - h;
422: geom.dimEmbed = dimEmbed;
423: for (f = 0, v = 0; f < Nf; ++f) {
424: void * const ctx = ctxs ? ctxs[f] : NULL;
426: if (!sp[f]) continue;
427: PetscDualSpaceGetDimension(sp[f], &spDim);
428: for (d = 0; d < spDim; ++d) {
429: if (funcs[f]) {
430: PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);
431: } else {
432: for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
433: }
434: v += numComp[f];
435: }
436: }
437: DMPlexVecSetClosure(dm, section, localX, p, values, mode);
438: }
439: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
440: if (h || !maxHeight) {
441: for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&sp[f]);}
442: }
443: }
444: PetscFree2(sp, numComp);
445: if (maxHeight > 0) {
446: for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&cellsp[f]);}
447: PetscFree(cellsp);
448: }
449: return(0);
450: }
454: /*@C
455: DMPlexProjectFunction - This projects the given function into the function space provided.
457: Input Parameters:
458: + dm - The DM
459: . funcs - The coordinate functions to evaluate, one per field
460: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
461: - mode - The insertion mode for values
463: Output Parameter:
464: . X - vector
466: Calling sequence of func:
467: $ func(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);
469: + dim - The spatial dimension
470: . x - The coordinates
471: . Nf - The number of fields
472: . u - The output field values
473: - ctx - optional user-defined function context
475: Level: developer
477: .seealso: DMPlexComputeL2Diff()
478: @*/
479: PetscErrorCode DMPlexProjectFunction(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
480: {
481: Vec localX;
486: DMGetLocalVector(dm, &localX);
487: DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
488: DMLocalToGlobalBegin(dm, localX, mode, X);
489: DMLocalToGlobalEnd(dm, localX, mode, X);
490: DMRestoreLocalVector(dm, &localX);
491: return(0);
492: }
496: PetscErrorCode DMPlexProjectFieldLocal(DM dm, Vec localU, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec localX)
497: {
498: DM dmAux;
499: PetscDS prob, probAux = NULL;
500: Vec A;
501: PetscSection section, sectionAux = NULL;
502: PetscDualSpace *sp;
503: PetscInt *Ncf;
504: PetscScalar *values, *u, *u_x, *a, *a_x;
505: PetscReal *x, *v0, *J, *invJ, detJ;
506: PetscInt Nf, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
507: PetscErrorCode ierr;
510: DMPlexGetMaxProjectionHeight(dm,&maxHeight);
511: if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");}
512: DMGetDS(dm, &prob);
513: DMGetDimension(dm, &dim);
514: DMGetDefaultSection(dm, §ion);
515: PetscSectionGetNumFields(section, &Nf);
516: PetscMalloc2(Nf, &sp, Nf, &Ncf);
517: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
518: PetscDSGetTotalDimension(prob, &totDim);
519: PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
520: PetscDSGetRefCoordArrays(prob, &x, NULL);
521: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
522: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
523: if (dmAux) {
524: DMGetDS(dmAux, &probAux);
525: DMGetDefaultSection(dmAux, §ionAux);
526: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
527: }
528: DMPlexInsertBoundaryValues(dm, localU, 0.0, NULL, NULL, NULL);
529: DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
530: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
531: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
532: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
533: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
534: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
535: for (c = cStart; c < cEnd; ++c) {
536: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
538: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
539: DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
540: if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
541: for (f = 0, v = 0; f < Nf; ++f) {
542: PetscObject obj;
543: PetscClassId id;
545: PetscDSGetDiscretization(prob, f, &obj);
546: PetscObjectGetClassId(obj, &id);
547: if (id == PETSCFE_CLASSID) {
548: PetscFE fe = (PetscFE) obj;
550: PetscFEGetDualSpace(fe, &sp[f]);
551: PetscObjectReference((PetscObject) sp[f]);
552: PetscFEGetNumComponents(fe, &Ncf[f]);
553: } else if (id == PETSCFV_CLASSID) {
554: PetscFV fv = (PetscFV) obj;
556: PetscFVGetNumComponents(fv, &Ncf[f]);
557: PetscFVGetDualSpace(fv, &sp[f]);
558: PetscObjectReference((PetscObject) sp[f]);
559: }
560: PetscDualSpaceGetDimension(sp[f], &spDim);
561: for (d = 0; d < spDim; ++d) {
562: PetscQuadrature quad;
563: const PetscReal *points, *weights;
564: PetscInt numPoints, q;
566: if (funcs[f]) {
567: PetscDualSpaceGetFunctional(sp[f], d, &quad);
568: PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
569: for (q = 0; q < numPoints; ++q) {
570: CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
571: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);
572: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);
573: (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]);
574: }
575: } else {
576: for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0;
577: }
578: v += Ncf[f];
579: }
580: }
581: DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);
582: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
583: DMPlexVecSetClosure(dm, section, localX, c, values, mode);
584: }
585: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
586: PetscFree3(v0,J,invJ);
587: for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&sp[f]);}
588: PetscFree2(sp, Ncf);
589: return(0);
590: }
594: static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX)
595: {
596: PetscErrorCode (**funcs)(PetscInt, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
597: void **ctxs;
598: PetscInt numFields;
599: PetscErrorCode ierr;
602: DMGetNumFields(dm, &numFields);
603: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
604: funcs[field] = func;
605: ctxs[field] = ctx;
606: DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);
607: PetscFree2(funcs,ctxs);
608: return(0);
609: }
613: /* This ignores numcomps/comps */
614: static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
615: PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
616: {
617: PetscDS prob;
618: PetscSF sf;
619: DM dmFace, dmCell, dmGrad;
620: const PetscScalar *facegeom, *cellgeom, *grad;
621: const PetscInt *leaves;
622: PetscScalar *x, *fx;
623: PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i;
624: PetscErrorCode ierr;
627: DMGetPointSF(dm, &sf);
628: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
629: nleaves = PetscMax(0, nleaves);
630: DMGetDimension(dm, &dim);
631: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
632: DMGetDS(dm, &prob);
633: VecGetDM(faceGeometry, &dmFace);
634: VecGetArrayRead(faceGeometry, &facegeom);
635: VecGetDM(cellGeometry, &dmCell);
636: VecGetArrayRead(cellGeometry, &cellgeom);
637: if (Grad) {
638: PetscFV fv;
640: PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
641: VecGetDM(Grad, &dmGrad);
642: VecGetArrayRead(Grad, &grad);
643: PetscFVGetNumComponents(fv, &pdim);
644: DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
645: }
646: VecGetArray(locX, &x);
647: for (i = 0; i < numids; ++i) {
648: IS faceIS;
649: const PetscInt *faces;
650: PetscInt numFaces, f;
652: DMLabelGetStratumIS(label, ids[i], &faceIS);
653: if (!faceIS) continue; /* No points with that id on this process */
654: ISGetLocalSize(faceIS, &numFaces);
655: ISGetIndices(faceIS, &faces);
656: for (f = 0; f < numFaces; ++f) {
657: const PetscInt face = faces[f], *cells;
658: const PetscFVFaceGeom *fg;
660: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
661: PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
662: if (loc >= 0) continue;
663: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
664: DMPlexGetSupport(dm, face, &cells);
665: if (Grad) {
666: const PetscFVCellGeom *cg;
667: const PetscScalar *cx, *cgrad;
668: PetscScalar *xG;
669: PetscReal dx[3];
670: PetscInt d;
672: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
673: DMPlexPointLocalRead(dm, cells[0], x, &cx);
674: DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
675: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
676: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
677: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
678: (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
679: } else {
680: const PetscScalar *xI;
681: PetscScalar *xG;
683: DMPlexPointLocalRead(dm, cells[0], x, &xI);
684: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
685: (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
686: }
687: }
688: ISRestoreIndices(faceIS, &faces);
689: ISDestroy(&faceIS);
690: }
691: VecRestoreArray(locX, &x);
692: if (Grad) {
693: DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
694: VecRestoreArrayRead(Grad, &grad);
695: }
696: VecRestoreArrayRead(cellGeometry, &cellgeom);
697: VecRestoreArrayRead(faceGeometry, &facegeom);
698: return(0);
699: }
703: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
704: {
705: PetscInt numBd, b;
714: DMPlexGetNumBoundary(dm, &numBd);
715: for (b = 0; b < numBd; ++b) {
716: PetscBool isEssential;
717: const char *labelname;
718: DMLabel label;
719: PetscInt field;
720: PetscObject obj;
721: PetscClassId id;
722: void (*func)();
723: PetscInt numids;
724: const PetscInt *ids;
725: void *ctx;
727: DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);
728: if (!isEssential) continue;
729: DMPlexGetLabel(dm, labelname, &label);
730: DMGetField(dm, field, &obj);
731: PetscObjectGetClassId(obj, &id);
732: if (id == PETSCFE_CLASSID) {
733: DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
734: } else if (id == PETSCFV_CLASSID) {
735: if (!faceGeomFVM) continue;
736: DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
737: field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
738: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
739: }
740: return(0);
741: }
745: /*@C
746: DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
748: Input Parameters:
749: + dm - The DM
750: . funcs - The functions to evaluate for each field component
751: . ctxs - Optional array of contexts to pass to each function, or NULL.
752: - X - The coefficient vector u_h
754: Output Parameter:
755: . diff - The diff ||u - u_h||_2
757: Level: developer
759: .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff()
760: @*/
761: PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
762: {
763: const PetscInt debug = 0;
764: PetscSection section;
765: PetscQuadrature quad;
766: Vec localX;
767: PetscScalar *funcVal, *interpolant;
768: PetscReal *coords, *v0, *J, *invJ, detJ;
769: PetscReal localDiff = 0.0;
770: const PetscReal *quadPoints, *quadWeights;
771: PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
772: PetscErrorCode ierr;
775: DMGetDimension(dm, &dim);
776: DMGetDefaultSection(dm, §ion);
777: PetscSectionGetNumFields(section, &numFields);
778: DMGetLocalVector(dm, &localX);
779: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
780: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
781: for (field = 0; field < numFields; ++field) {
782: PetscObject obj;
783: PetscClassId id;
784: PetscInt Nc;
786: DMGetField(dm, field, &obj);
787: PetscObjectGetClassId(obj, &id);
788: if (id == PETSCFE_CLASSID) {
789: PetscFE fe = (PetscFE) obj;
791: PetscFEGetQuadrature(fe, &quad);
792: PetscFEGetNumComponents(fe, &Nc);
793: } else if (id == PETSCFV_CLASSID) {
794: PetscFV fv = (PetscFV) obj;
796: PetscFVGetQuadrature(fv, &quad);
797: PetscFVGetNumComponents(fv, &Nc);
798: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
799: numComponents += Nc;
800: }
801: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
802: DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
803: PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
804: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
805: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
806: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
807: for (c = cStart; c < cEnd; ++c) {
808: PetscScalar *x = NULL;
809: PetscReal elemDiff = 0.0;
811: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
812: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
813: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
815: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
816: PetscObject obj;
817: PetscClassId id;
818: void * const ctx = ctxs ? ctxs[field] : NULL;
819: PetscInt Nb, Nc, q, fc;
821: DMGetField(dm, field, &obj);
822: PetscObjectGetClassId(obj, &id);
823: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
824: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
825: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
826: if (debug) {
827: char title[1024];
828: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
829: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
830: }
831: for (q = 0; q < numQuadPoints; ++q) {
832: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
833: (*funcs[field])(dim, coords, numFields, funcVal, ctx);
834: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
835: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
836: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
837: for (fc = 0; fc < Nc; ++fc) {
838: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
839: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
840: }
841: }
842: fieldOffset += Nb*Nc;
843: }
844: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
845: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
846: localDiff += elemDiff;
847: }
848: PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
849: DMRestoreLocalVector(dm, &localX);
850: MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
851: *diff = PetscSqrtReal(*diff);
852: return(0);
853: }
857: /*@C
858: DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
860: Input Parameters:
861: + dm - The DM
862: . funcs - The gradient functions to evaluate for each field component
863: . ctxs - Optional array of contexts to pass to each function, or NULL.
864: . X - The coefficient vector u_h
865: - n - The vector to project along
867: Output Parameter:
868: . diff - The diff ||(grad u - grad u_h) . n||_2
870: Level: developer
872: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
873: @*/
874: PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
875: {
876: const PetscInt debug = 0;
877: PetscSection section;
878: PetscQuadrature quad;
879: Vec localX;
880: PetscScalar *funcVal, *interpolantVec;
881: PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
882: PetscReal localDiff = 0.0;
883: PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
884: PetscErrorCode ierr;
887: DMGetDimension(dm, &dim);
888: DMGetDefaultSection(dm, §ion);
889: PetscSectionGetNumFields(section, &numFields);
890: DMGetLocalVector(dm, &localX);
891: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
892: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
893: for (field = 0; field < numFields; ++field) {
894: PetscFE fe;
895: PetscInt Nc;
897: DMGetField(dm, field, (PetscObject *) &fe);
898: PetscFEGetQuadrature(fe, &quad);
899: PetscFEGetNumComponents(fe, &Nc);
900: numComponents += Nc;
901: }
902: /* DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
903: PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
904: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
905: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
906: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
907: for (c = cStart; c < cEnd; ++c) {
908: PetscScalar *x = NULL;
909: PetscReal elemDiff = 0.0;
911: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
912: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
913: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
915: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
916: PetscFE fe;
917: void * const ctx = ctxs ? ctxs[field] : NULL;
918: const PetscReal *quadPoints, *quadWeights;
919: PetscReal *basisDer;
920: PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
922: DMGetField(dm, field, (PetscObject *) &fe);
923: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
924: PetscFEGetDimension(fe, &Nb);
925: PetscFEGetNumComponents(fe, &Ncomp);
926: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
927: if (debug) {
928: char title[1024];
929: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
930: DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);
931: }
932: for (q = 0; q < numQuadPoints; ++q) {
933: for (d = 0; d < dim; d++) {
934: coords[d] = v0[d];
935: for (e = 0; e < dim; e++) {
936: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
937: }
938: }
939: (*funcs[field])(dim, coords, n, numFields, funcVal, ctx);
940: for (fc = 0; fc < Ncomp; ++fc) {
941: PetscScalar interpolant = 0.0;
943: for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
944: for (f = 0; f < Nb; ++f) {
945: const PetscInt fidx = f*Ncomp+fc;
947: for (d = 0; d < dim; ++d) {
948: realSpaceDer[d] = 0.0;
949: for (g = 0; g < dim; ++g) {
950: realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
951: }
952: interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
953: }
954: }
955: for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
956: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
957: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
958: }
959: }
960: comp += Ncomp;
961: fieldOffset += Nb*Ncomp;
962: }
963: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
964: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
965: localDiff += elemDiff;
966: }
967: PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
968: DMRestoreLocalVector(dm, &localX);
969: MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
970: *diff = PetscSqrtReal(*diff);
971: return(0);
972: }
976: /*@C
977: DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
979: Input Parameters:
980: + dm - The DM
981: . funcs - The functions to evaluate for each field component
982: . ctxs - Optional array of contexts to pass to each function, or NULL.
983: - X - The coefficient vector u_h
985: Output Parameter:
986: . diff - The array of differences, ||u^f - u^f_h||_2
988: Level: developer
990: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff()
991: @*/
992: PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
993: {
994: const PetscInt debug = 0;
995: PetscSection section;
996: PetscQuadrature quad;
997: Vec localX;
998: PetscScalar *funcVal, *interpolant;
999: PetscReal *coords, *v0, *J, *invJ, detJ;
1000: PetscReal *localDiff;
1001: const PetscReal *quadPoints, *quadWeights;
1002: PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1003: PetscErrorCode ierr;
1006: DMGetDimension(dm, &dim);
1007: DMGetDefaultSection(dm, §ion);
1008: PetscSectionGetNumFields(section, &numFields);
1009: DMGetLocalVector(dm, &localX);
1010: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1011: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1012: for (field = 0; field < numFields; ++field) {
1013: PetscObject obj;
1014: PetscClassId id;
1015: PetscInt Nc;
1017: DMGetField(dm, field, &obj);
1018: PetscObjectGetClassId(obj, &id);
1019: if (id == PETSCFE_CLASSID) {
1020: PetscFE fe = (PetscFE) obj;
1022: PetscFEGetQuadrature(fe, &quad);
1023: PetscFEGetNumComponents(fe, &Nc);
1024: } else if (id == PETSCFV_CLASSID) {
1025: PetscFV fv = (PetscFV) obj;
1027: PetscFVGetQuadrature(fv, &quad);
1028: PetscFVGetNumComponents(fv, &Nc);
1029: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1030: numComponents += Nc;
1031: }
1032: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1033: DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
1034: PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1035: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1036: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1037: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1038: for (c = cStart; c < cEnd; ++c) {
1039: PetscScalar *x = NULL;
1041: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
1042: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1043: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1045: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1046: PetscObject obj;
1047: PetscClassId id;
1048: void * const ctx = ctxs ? ctxs[field] : NULL;
1049: PetscInt Nb, Nc, q, fc;
1051: PetscReal elemDiff = 0.0;
1053: DMGetField(dm, field, &obj);
1054: PetscObjectGetClassId(obj, &id);
1055: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1056: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1057: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1058: if (debug) {
1059: char title[1024];
1060: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1061: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1062: }
1063: for (q = 0; q < numQuadPoints; ++q) {
1064: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1065: (*funcs[field])(dim, coords, numFields, funcVal, ctx);
1066: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1067: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1068: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1069: for (fc = 0; fc < Nc; ++fc) {
1070: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
1071: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1072: }
1073: }
1074: fieldOffset += Nb*Nc;
1075: localDiff[field] += elemDiff;
1076: }
1077: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1078: }
1079: DMRestoreLocalVector(dm, &localX);
1080: MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1081: for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1082: PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);
1083: return(0);
1084: }
1088: /*@
1089: DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user
1091: Input Parameters:
1092: + dm - The mesh
1093: . X - Local input vector
1094: - user - The user context
1096: Output Parameter:
1097: . integral - Local integral for each field
1099: Level: developer
1101: .seealso: DMPlexComputeResidualFEM()
1102: @*/
1103: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1104: {
1105: DM_Plex *mesh = (DM_Plex *) dm->data;
1106: DM dmAux;
1107: Vec localX, A;
1108: PetscDS prob, probAux = NULL;
1109: PetscSection section, sectionAux;
1110: PetscFECellGeom *cgeom;
1111: PetscScalar *u, *a = NULL;
1112: PetscReal *lintegral, *vol;
1113: PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c;
1114: PetscInt totDim, totDimAux;
1115: PetscErrorCode ierr;
1118: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1119: DMGetLocalVector(dm, &localX);
1120: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1121: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1122: DMGetDimension(dm, &dim);
1123: DMGetDefaultSection(dm, §ion);
1124: DMGetDS(dm, &prob);
1125: PetscDSGetTotalDimension(prob, &totDim);
1126: PetscSectionGetNumFields(section, &Nf);
1127: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1128: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1129: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1130: numCells = cEnd - cStart;
1131: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1132: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1133: if (dmAux) {
1134: DMGetDefaultSection(dmAux, §ionAux);
1135: DMGetDS(dmAux, &probAux);
1136: PetscDSGetTotalDimension(probAux, &totDimAux);
1137: }
1138: DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);
1139: PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);
1140: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1141: for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1142: for (c = cStart; c < cEnd; ++c) {
1143: PetscScalar *x = NULL;
1144: PetscInt i;
1146: DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);
1147: DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);
1148: if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1149: DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1150: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1151: DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1152: if (dmAux) {
1153: DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1154: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1155: DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1156: }
1157: }
1158: for (f = 0; f < Nf; ++f) {
1159: PetscObject obj;
1160: PetscClassId id;
1161: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1163: PetscDSGetDiscretization(prob, f, &obj);
1164: PetscObjectGetClassId(obj, &id);
1165: if (id == PETSCFE_CLASSID) {
1166: PetscFE fe = (PetscFE) obj;
1167: PetscQuadrature q;
1168: PetscInt Nq, Nb;
1170: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1171: PetscFEGetQuadrature(fe, &q);
1172: PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1173: PetscFEGetDimension(fe, &Nb);
1174: blockSize = Nb*Nq;
1175: batchSize = numBlocks * blockSize;
1176: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1177: numChunks = numCells / (numBatches*batchSize);
1178: Ne = numChunks*numBatches*batchSize;
1179: Nr = numCells % (numBatches*batchSize);
1180: offset = numCells - Nr;
1181: PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);
1182: PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1183: } else if (id == PETSCFV_CLASSID) {
1184: /* PetscFV fv = (PetscFV) obj; */
1185: PetscInt foff;
1186: PetscPointFunc obj_func;
1187: PetscScalar lint;
1189: PetscDSGetObjective(prob, f, &obj_func);
1190: PetscDSGetFieldOffset(prob, f, &foff);
1191: if (obj_func) {
1192: for (c = 0; c < numCells; ++c) {
1193: /* TODO: Need full pointwise interpolation and get centroid for x argument */
1194: obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint);
1195: lintegral[f] = PetscRealPart(lint)*vol[c];
1196: }
1197: }
1198: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1199: }
1200: if (dmAux) {PetscFree(a);}
1201: if (mesh->printFEM) {
1202: PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1203: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1204: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1205: }
1206: DMRestoreLocalVector(dm, &localX);
1207: MPI_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1208: PetscFree4(lintegral,u,cgeom,vol);
1209: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1210: return(0);
1211: }
1215: /*@
1216: DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1218: Input Parameters:
1219: + dmf - The fine mesh
1220: . dmc - The coarse mesh
1221: - user - The user context
1223: Output Parameter:
1224: . In - The interpolation matrix
1226: Note:
1227: The first member of the user context must be an FEMContext.
1229: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1230: like a GPU, or vectorize on a multicore machine.
1232: Level: developer
1234: .seealso: DMPlexComputeJacobianFEM()
1235: @*/
1236: PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1237: {
1238: DM_Plex *mesh = (DM_Plex *) dmc->data;
1239: const char *name = "Interpolator";
1240: PetscDS prob;
1241: PetscFE *feRef;
1242: PetscFV *fvRef;
1243: PetscSection fsection, fglobalSection;
1244: PetscSection csection, cglobalSection;
1245: PetscScalar *elemMat;
1246: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1247: PetscInt cTotDim, rTotDim = 0;
1248: PetscErrorCode ierr;
1251: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1252: DMGetDimension(dmf, &dim);
1253: DMGetDefaultSection(dmf, &fsection);
1254: DMGetDefaultGlobalSection(dmf, &fglobalSection);
1255: DMGetDefaultSection(dmc, &csection);
1256: DMGetDefaultGlobalSection(dmc, &cglobalSection);
1257: PetscSectionGetNumFields(fsection, &Nf);
1258: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1259: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1260: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1261: DMGetDS(dmf, &prob);
1262: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1263: for (f = 0; f < Nf; ++f) {
1264: PetscObject obj;
1265: PetscClassId id;
1266: PetscInt rNb = 0, Nc = 0;
1268: PetscDSGetDiscretization(prob, f, &obj);
1269: PetscObjectGetClassId(obj, &id);
1270: if (id == PETSCFE_CLASSID) {
1271: PetscFE fe = (PetscFE) obj;
1273: PetscFERefine(fe, &feRef[f]);
1274: PetscFEGetDimension(feRef[f], &rNb);
1275: PetscFEGetNumComponents(fe, &Nc);
1276: } else if (id == PETSCFV_CLASSID) {
1277: PetscFV fv = (PetscFV) obj;
1278: PetscDualSpace Q;
1280: PetscFVRefine(fv, &fvRef[f]);
1281: PetscFVGetDualSpace(fvRef[f], &Q);
1282: PetscDualSpaceGetDimension(Q, &rNb);
1283: PetscFVGetNumComponents(fv, &Nc);
1284: }
1285: rTotDim += rNb*Nc;
1286: }
1287: PetscDSGetTotalDimension(prob, &cTotDim);
1288: PetscMalloc1(rTotDim*cTotDim,&elemMat);
1289: PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1290: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1291: PetscDualSpace Qref;
1292: PetscQuadrature f;
1293: const PetscReal *qpoints, *qweights;
1294: PetscReal *points;
1295: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
1297: /* Compose points from all dual basis functionals */
1298: if (feRef[fieldI]) {
1299: PetscFEGetDualSpace(feRef[fieldI], &Qref);
1300: PetscFEGetNumComponents(feRef[fieldI], &Nc);
1301: } else {
1302: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1303: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1304: }
1305: PetscDualSpaceGetDimension(Qref, &fpdim);
1306: for (i = 0; i < fpdim; ++i) {
1307: PetscDualSpaceGetFunctional(Qref, i, &f);
1308: PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1309: npoints += Np;
1310: }
1311: PetscMalloc1(npoints*dim,&points);
1312: for (i = 0, k = 0; i < fpdim; ++i) {
1313: PetscDualSpaceGetFunctional(Qref, i, &f);
1314: PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1315: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1316: }
1318: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1319: PetscObject obj;
1320: PetscClassId id;
1321: PetscReal *B;
1322: PetscInt NcJ = 0, cpdim = 0, j;
1324: PetscDSGetDiscretization(prob, fieldJ, &obj);
1325: PetscObjectGetClassId(obj, &id);
1326: if (id == PETSCFE_CLASSID) {
1327: PetscFE fe = (PetscFE) obj;
1329: /* Evaluate basis at points */
1330: PetscFEGetNumComponents(fe, &NcJ);
1331: PetscFEGetDimension(fe, &cpdim);
1332: /* For now, fields only interpolate themselves */
1333: if (fieldI == fieldJ) {
1334: if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ);
1335: PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1336: for (i = 0, k = 0; i < fpdim; ++i) {
1337: PetscDualSpaceGetFunctional(Qref, i, &f);
1338: PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1339: for (p = 0; p < Np; ++p, ++k) {
1340: for (j = 0; j < cpdim; ++j) {
1341: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p];
1342: }
1343: }
1344: }
1345: PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1346: }
1347: } else if (id == PETSCFV_CLASSID) {
1348: PetscFV fv = (PetscFV) obj;
1350: /* Evaluate constant function at points */
1351: PetscFVGetNumComponents(fv, &NcJ);
1352: cpdim = 1;
1353: /* For now, fields only interpolate themselves */
1354: if (fieldI == fieldJ) {
1355: if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ);
1356: for (i = 0, k = 0; i < fpdim; ++i) {
1357: PetscDualSpaceGetFunctional(Qref, i, &f);
1358: PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1359: for (p = 0; p < Np; ++p, ++k) {
1360: for (j = 0; j < cpdim; ++j) {
1361: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1362: }
1363: }
1364: }
1365: }
1366: }
1367: offsetJ += cpdim*NcJ;
1368: }
1369: offsetI += fpdim*Nc;
1370: PetscFree(points);
1371: }
1372: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1373: /* Preallocate matrix */
1374: {
1375: PetscHashJK ht;
1376: PetscLayout rLayout;
1377: PetscInt *dnz, *onz, *cellCIndices, *cellFIndices;
1378: PetscInt locRows, rStart, rEnd, cell, r;
1380: MatGetLocalSize(In, &locRows, NULL);
1381: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1382: PetscLayoutSetLocalSize(rLayout, locRows);
1383: PetscLayoutSetBlockSize(rLayout, 1);
1384: PetscLayoutSetUp(rLayout);
1385: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1386: PetscLayoutDestroy(&rLayout);
1387: PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1388: PetscHashJKCreate(&ht);
1389: for (cell = cStart; cell < cEnd; ++cell) {
1390: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1391: for (r = 0; r < rTotDim; ++r) {
1392: PetscHashJKKey key;
1393: PetscHashJKIter missing, iter;
1395: key.j = cellFIndices[r];
1396: if (key.j < 0) continue;
1397: for (c = 0; c < cTotDim; ++c) {
1398: key.k = cellCIndices[c];
1399: if (key.k < 0) continue;
1400: PetscHashJKPut(ht, key, &missing, &iter);
1401: if (missing) {
1402: PetscHashJKSet(ht, iter, 1);
1403: if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1404: else ++onz[key.j-rStart];
1405: }
1406: }
1407: }
1408: }
1409: PetscHashJKDestroy(&ht);
1410: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1411: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1412: PetscFree4(dnz,onz,cellCIndices,cellFIndices);
1413: }
1414: /* Fill matrix */
1415: MatZeroEntries(In);
1416: for (c = cStart; c < cEnd; ++c) {
1417: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1418: }
1419: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1420: PetscFree2(feRef,fvRef);
1421: PetscFree(elemMat);
1422: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1423: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1424: if (mesh->printFEM) {
1425: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1426: MatChop(In, 1.0e-10);
1427: MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1428: }
1429: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1430: return(0);
1431: }
1435: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1436: {
1437: PetscDS prob;
1438: PetscFE *feRef;
1439: PetscFV *fvRef;
1440: Vec fv, cv;
1441: IS fis, cis;
1442: PetscSection fsection, fglobalSection, csection, cglobalSection;
1443: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1444: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;
1448: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1449: DMGetDimension(dmf, &dim);
1450: DMGetDefaultSection(dmf, &fsection);
1451: DMGetDefaultGlobalSection(dmf, &fglobalSection);
1452: DMGetDefaultSection(dmc, &csection);
1453: DMGetDefaultGlobalSection(dmc, &cglobalSection);
1454: PetscSectionGetNumFields(fsection, &Nf);
1455: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1456: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1457: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1458: DMGetDS(dmc, &prob);
1459: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1460: for (f = 0; f < Nf; ++f) {
1461: PetscObject obj;
1462: PetscClassId id;
1463: PetscInt fNb = 0, Nc = 0;
1465: PetscDSGetDiscretization(prob, f, &obj);
1466: PetscObjectGetClassId(obj, &id);
1467: if (id == PETSCFE_CLASSID) {
1468: PetscFE fe = (PetscFE) obj;
1470: PetscFERefine(fe, &feRef[f]);
1471: PetscFEGetDimension(feRef[f], &fNb);
1472: PetscFEGetNumComponents(fe, &Nc);
1473: } else if (id == PETSCFV_CLASSID) {
1474: PetscFV fv = (PetscFV) obj;
1475: PetscDualSpace Q;
1477: PetscFVRefine(fv, &fvRef[f]);
1478: PetscFVGetDualSpace(fvRef[f], &Q);
1479: PetscDualSpaceGetDimension(Q, &fNb);
1480: PetscFVGetNumComponents(fv, &Nc);
1481: }
1482: fTotDim += fNb*Nc;
1483: }
1484: PetscDSGetTotalDimension(prob, &cTotDim);
1485: PetscMalloc1(cTotDim,&cmap);
1486: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1487: PetscFE feC;
1488: PetscFV fvC;
1489: PetscDualSpace QF, QC;
1490: PetscInt NcF, NcC, fpdim, cpdim;
1492: if (feRef[field]) {
1493: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1494: PetscFEGetNumComponents(feC, &NcC);
1495: PetscFEGetNumComponents(feRef[field], &NcF);
1496: PetscFEGetDualSpace(feRef[field], &QF);
1497: PetscDualSpaceGetDimension(QF, &fpdim);
1498: PetscFEGetDualSpace(feC, &QC);
1499: PetscDualSpaceGetDimension(QC, &cpdim);
1500: } else {
1501: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1502: PetscFVGetNumComponents(fvC, &NcC);
1503: PetscFVGetNumComponents(fvRef[field], &NcF);
1504: PetscFVGetDualSpace(fvRef[field], &QF);
1505: PetscDualSpaceGetDimension(QF, &fpdim);
1506: PetscFVGetDualSpace(fvC, &QC);
1507: PetscDualSpaceGetDimension(QC, &cpdim);
1508: }
1509: if (NcF != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", NcF, NcC);
1510: for (c = 0; c < cpdim; ++c) {
1511: PetscQuadrature cfunc;
1512: const PetscReal *cqpoints;
1513: PetscInt NpC;
1514: PetscBool found = PETSC_FALSE;
1516: PetscDualSpaceGetFunctional(QC, c, &cfunc);
1517: PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1518: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1519: for (f = 0; f < fpdim; ++f) {
1520: PetscQuadrature ffunc;
1521: const PetscReal *fqpoints;
1522: PetscReal sum = 0.0;
1523: PetscInt NpF, comp;
1525: PetscDualSpaceGetFunctional(QF, f, &ffunc);
1526: PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1527: if (NpC != NpF) continue;
1528: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1529: if (sum > 1.0e-9) continue;
1530: for (comp = 0; comp < NcC; ++comp) {
1531: cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1532: }
1533: found = PETSC_TRUE;
1534: break;
1535: }
1536: if (!found) {
1537: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1538: if (fvRef[field]) {
1539: PetscInt comp;
1540: for (comp = 0; comp < NcC; ++comp) {
1541: cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1542: }
1543: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1544: }
1545: }
1546: offsetC += cpdim*NcC;
1547: offsetF += fpdim*NcF;
1548: }
1549: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1550: PetscFree2(feRef,fvRef);
1552: DMGetGlobalVector(dmf, &fv);
1553: DMGetGlobalVector(dmc, &cv);
1554: VecGetOwnershipRange(cv, &startC, NULL);
1555: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1556: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1557: PetscMalloc1(m,&cindices);
1558: PetscMalloc1(m,&findices);
1559: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1560: for (c = cStart; c < cEnd; ++c) {
1561: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1562: for (d = 0; d < cTotDim; ++d) {
1563: if (cellCIndices[d] < 0) continue;
1564: if ((findices[cellCIndices[d]-startC] >= 0) && (findices[cellCIndices[d]-startC] != cellFIndices[cmap[d]])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coarse dof %d maps to both %d and %d", cindices[cellCIndices[d]-startC], findices[cellCIndices[d]-startC], cellFIndices[cmap[d]]);
1565: cindices[cellCIndices[d]-startC] = cellCIndices[d];
1566: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1567: }
1568: }
1569: PetscFree(cmap);
1570: PetscFree2(cellCIndices,cellFIndices);
1572: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1573: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1574: VecScatterCreate(cv, cis, fv, fis, sc);
1575: ISDestroy(&cis);
1576: ISDestroy(&fis);
1577: DMRestoreGlobalVector(dmf, &fv);
1578: DMRestoreGlobalVector(dmc, &cv);
1579: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1580: return(0);
1581: }