Actual source code: plexfem.c
petsc-3.7.5 2017-01-01
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, PetscReal t, 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, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody;
159: ctx[0] = dimEmbed;
160: ctx[1] = d;
161: DMProjectFunction(dm, 0.0, &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(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
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(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
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 DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, 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: } else {
278: PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
279: if (!sp[f]) continue;
280: }
281: } else if (id == PETSCFV_CLASSID) {
282: PetscFV fv = (PetscFV) obj;
284: PetscFVGetNumComponents(fv, &numComp[f]);
285: PetscFVGetDualSpace(fv, &sp[f]);
286: PetscObjectReference((PetscObject) sp[f]);
287: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
288: PetscDualSpaceGetDimension(sp[f], &spDim);
289: totDim += spDim*numComp[f];
290: }
291: DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
292: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
293: if (!totDim) continue;
294: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
295: DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
296: for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
297: for (i = 0; i < numIds; ++i) {
298: IS pointIS;
299: const PetscInt *points;
300: PetscInt n, p;
302: DMLabelGetStratumIS(label, ids[i], &pointIS);
303: if (!pointIS) continue; /* No points with that id on this process */
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, time, &geom, numComp[f], funcs[f], ctx, &values[v]);
322: if (ierr) {
323: PetscErrorCode ierr2;
324: ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2);
325: ierr2 = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2);
326:
327: }
328: } else {
329: for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
330: }
331: v += numComp[f];
332: }
333: }
334: DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);
335: }
336: ISRestoreIndices(pointIS, &points);
337: ISDestroy(&pointIS);
338: }
339: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
340: DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
341: }
342: PetscFree2(sp, numComp);
343: if (maxHeight > 0) {
344: PetscFree(cellsp);
345: }
346: return(0);
347: }
351: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
352: {
353: PetscDualSpace *sp, *cellsp;
354: PetscInt *numComp;
355: PetscSection section;
356: PetscScalar *values;
357: PetscInt Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight;
358: PetscBool *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE;
359: PetscErrorCode ierr;
362: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
363: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
364: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
365: DMGetDefaultSection(dm, §ion);
366: PetscSectionGetNumFields(section, &Nf);
367: PetscMalloc3(Nf, &isFE, Nf, &sp, Nf, &numComp);
368: DMGetDimension(dm, &dim);
369: DMGetCoordinateDim(dm, &dimEmbed);
370: DMPlexGetMaxProjectionHeight(dm,&maxHeight);
371: if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
372: if (maxHeight > 0) {
373: PetscMalloc1(Nf,&cellsp);
374: }
375: else {
376: cellsp = sp;
377: }
378: for (f = 0; f < Nf; ++f) {
379: PetscObject obj;
380: PetscClassId id;
382: DMGetField(dm, f, &obj);
383: PetscObjectGetClassId(obj, &id);
384: if (id == PETSCFE_CLASSID) {
385: PetscFE fe = (PetscFE) obj;
387: hasFE = PETSC_TRUE;
388: isFE[f] = PETSC_TRUE;
389: PetscFEGetNumComponents(fe, &numComp[f]);
390: PetscFEGetDualSpace(fe, &cellsp[f]);
391: } else if (id == PETSCFV_CLASSID) {
392: PetscFV fv = (PetscFV) obj;
394: hasFV = PETSC_TRUE;
395: isFE[f] = PETSC_FALSE;
396: PetscFVGetNumComponents(fv, &numComp[f]);
397: PetscFVGetDualSpace(fv, &cellsp[f]);
398: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
399: }
400: for (h = 0; h <= maxHeight; h++) {
401: DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
402: if (!h) {pStart = cStart; pEnd = cEnd;}
403: if (pEnd <= pStart) continue;
404: totDim = 0;
405: for (f = 0; f < Nf; ++f) {
406: if (!h) {
407: sp[f] = cellsp[f];
408: }
409: else {
410: PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
411: if (!sp[f]) {
412: continue;
413: }
414: }
415: PetscDualSpaceGetDimension(sp[f], &spDim);
416: totDim += spDim*numComp[f];
417: }
418: DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
419: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
420: if (!totDim) continue;
421: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
422: for (p = pStart; p < pEnd; ++p) {
423: PetscFECellGeom fegeom;
424: PetscFVCellGeom fvgeom;
426: if (hasFE) {
427: DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);
428: fegeom.dim = dim - h;
429: fegeom.dimEmbed = dimEmbed;
430: }
431: if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
432: for (f = 0, v = 0; f < Nf; ++f) {
433: void * const ctx = ctxs ? ctxs[f] : NULL;
435: if (!sp[f]) continue;
436: PetscDualSpaceGetDimension(sp[f], &spDim);
437: for (d = 0; d < spDim; ++d) {
438: if (funcs[f]) {
439: if (isFE[f]) {PetscDualSpaceApply(sp[f], d, time, &fegeom, numComp[f], funcs[f], ctx, &values[v]);}
440: else {PetscDualSpaceApplyFVM(sp[f], d, time, &fvgeom, numComp[f], funcs[f], ctx, &values[v]);}
441: if (ierr) {
442: PetscErrorCode ierr2;
443: ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2);
444:
445: }
446: } else {
447: for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
448: }
449: v += numComp[f];
450: }
451: }
452: DMPlexVecSetClosure(dm, section, localX, p, values, mode);
453: }
454: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
455: }
456: PetscFree3(isFE, sp, numComp);
457: if (maxHeight > 0) {
458: PetscFree(cellsp);
459: }
460: return(0);
461: }
465: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, Vec localU,
466: void (**funcs)(PetscInt, PetscInt, PetscInt,
467: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
468: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
469: PetscReal, const PetscReal[], PetscScalar[]),
470: InsertMode mode, Vec localX)
471: {
472: DM dmAux;
473: PetscDS prob, probAux = NULL;
474: Vec A;
475: PetscSection section, sectionAux = NULL;
476: PetscDualSpace *sp;
477: PetscInt *Ncf;
478: PetscScalar *values, *u, *u_x, *a, *a_x;
479: PetscReal *x, *v0, *J, *invJ, detJ;
480: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
481: PetscInt Nf, NfAux = 0, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
482: PetscErrorCode ierr;
485: DMPlexGetMaxProjectionHeight(dm,&maxHeight);
486: if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");}
487: DMGetDS(dm, &prob);
488: DMGetDimension(dm, &dim);
489: DMGetDefaultSection(dm, §ion);
490: PetscSectionGetNumFields(section, &Nf);
491: PetscMalloc2(Nf, &sp, Nf, &Ncf);
492: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
493: PetscDSGetTotalDimension(prob, &totDim);
494: PetscDSGetComponentOffsets(prob, &uOff);
495: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
496: PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
497: PetscDSGetRefCoordArrays(prob, &x, NULL);
498: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
499: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
500: if (dmAux) {
501: DMGetDS(dmAux, &probAux);
502: PetscDSGetNumFields(probAux, &NfAux);
503: DMGetDefaultSection(dmAux, §ionAux);
504: PetscDSGetComponentOffsets(probAux, &aOff);
505: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
506: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
507: }
508: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, 0.0, NULL, NULL, NULL);
509: DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
510: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
511: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
512: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
513: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
514: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
515: for (c = cStart; c < cEnd; ++c) {
516: PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
518: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
519: DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
520: if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
521: for (f = 0, v = 0; f < Nf; ++f) {
522: PetscObject obj;
523: PetscClassId id;
525: PetscDSGetDiscretization(prob, f, &obj);
526: PetscObjectGetClassId(obj, &id);
527: if (id == PETSCFE_CLASSID) {
528: PetscFE fe = (PetscFE) obj;
530: PetscFEGetDualSpace(fe, &sp[f]);
531: PetscFEGetNumComponents(fe, &Ncf[f]);
532: } else if (id == PETSCFV_CLASSID) {
533: PetscFV fv = (PetscFV) obj;
535: PetscFVGetNumComponents(fv, &Ncf[f]);
536: PetscObjectReference((PetscObject) sp[f]);
537: }
538: PetscDualSpaceGetDimension(sp[f], &spDim);
539: for (d = 0; d < spDim; ++d) {
540: PetscQuadrature quad;
541: const PetscReal *points, *weights;
542: PetscInt numPoints, q;
544: if (funcs[f]) {
545: PetscDualSpaceGetFunctional(sp[f], d, &quad);
546: PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
547: for (q = 0; q < numPoints; ++q) {
548: CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
549: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);
550: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);
551: (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &values[v]);
552: }
553: } else {
554: for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0;
555: }
556: v += Ncf[f];
557: }
558: }
559: DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);
560: if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
561: DMPlexVecSetClosure(dm, section, localX, c, values, mode);
562: }
563: PetscFree3(v0,J,invJ);
564: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
565: PetscFree2(sp, Ncf);
566: return(0);
567: }
571: static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX)
572: {
573: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
574: void **ctxs;
575: PetscInt numFields;
576: PetscErrorCode ierr;
579: DMGetNumFields(dm, &numFields);
580: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
581: funcs[field] = func;
582: ctxs[field] = ctx;
583: DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);
584: PetscFree2(funcs,ctxs);
585: return(0);
586: }
590: /* This ignores numcomps/comps */
591: static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
592: PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
593: {
594: PetscDS prob;
595: PetscSF sf;
596: DM dmFace, dmCell, dmGrad;
597: const PetscScalar *facegeom, *cellgeom = NULL, *grad;
598: const PetscInt *leaves;
599: PetscScalar *x, *fx;
600: PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i;
601: PetscErrorCode ierr, ierru = 0;
604: DMGetPointSF(dm, &sf);
605: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
606: nleaves = PetscMax(0, nleaves);
607: DMGetDimension(dm, &dim);
608: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
609: DMGetDS(dm, &prob);
610: VecGetDM(faceGeometry, &dmFace);
611: VecGetArrayRead(faceGeometry, &facegeom);
612: if (cellGeometry) {
613: VecGetDM(cellGeometry, &dmCell);
614: VecGetArrayRead(cellGeometry, &cellgeom);
615: }
616: if (Grad) {
617: PetscFV fv;
619: PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
620: VecGetDM(Grad, &dmGrad);
621: VecGetArrayRead(Grad, &grad);
622: PetscFVGetNumComponents(fv, &pdim);
623: DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
624: }
625: VecGetArray(locX, &x);
626: for (i = 0; i < numids; ++i) {
627: IS faceIS;
628: const PetscInt *faces;
629: PetscInt numFaces, f;
631: DMLabelGetStratumIS(label, ids[i], &faceIS);
632: if (!faceIS) continue; /* No points with that id on this process */
633: ISGetLocalSize(faceIS, &numFaces);
634: ISGetIndices(faceIS, &faces);
635: for (f = 0; f < numFaces; ++f) {
636: const PetscInt face = faces[f], *cells;
637: PetscFVFaceGeom *fg;
639: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
640: PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
641: if (loc >= 0) continue;
642: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
643: DMPlexGetSupport(dm, face, &cells);
644: if (Grad) {
645: PetscFVCellGeom *cg;
646: PetscScalar *cx, *cgrad;
647: PetscScalar *xG;
648: PetscReal dx[3];
649: PetscInt d;
651: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
652: DMPlexPointLocalRead(dm, cells[0], x, &cx);
653: DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
654: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
655: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
656: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
657: ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
658: if (ierru) {
659: ISRestoreIndices(faceIS, &faces);
660: ISDestroy(&faceIS);
661: goto cleanup;
662: }
663: } else {
664: PetscScalar *xI;
665: PetscScalar *xG;
667: DMPlexPointLocalRead(dm, cells[0], x, &xI);
668: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
669: ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
670: if (ierru) {
671: ISRestoreIndices(faceIS, &faces);
672: ISDestroy(&faceIS);
673: goto cleanup;
674: }
675: }
676: }
677: ISRestoreIndices(faceIS, &faces);
678: ISDestroy(&faceIS);
679: }
680: cleanup:
681: VecRestoreArray(locX, &x);
682: if (Grad) {
683: DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
684: VecRestoreArrayRead(Grad, &grad);
685: }
686: if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
687: VecRestoreArrayRead(faceGeometry, &facegeom);
688: CHKERRQ(ierru);
689: return(0);
690: }
694: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
695: {
696: PetscInt numBd, b;
705: DMGetNumBoundary(dm, &numBd);
706: for (b = 0; b < numBd; ++b) {
707: PetscBool isEssential;
708: const char *labelname;
709: DMLabel label;
710: PetscInt field;
711: PetscObject obj;
712: PetscClassId id;
713: void (*func)();
714: PetscInt numids;
715: const PetscInt *ids;
716: void *ctx;
718: DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);
719: if (insertEssential != isEssential) continue;
720: DMGetLabel(dm, labelname, &label);
721: DMGetField(dm, field, &obj);
722: PetscObjectGetClassId(obj, &id);
723: if (id == PETSCFE_CLASSID) {
724: if (!isEssential) continue; /* for FEM, there is no insertion to be done for non-essential boundary conditions */
725: DMPlexLabelAddCells(dm,label);
726: DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
727: DMPlexLabelClearCells(dm,label);
728: } else if (id == PETSCFV_CLASSID) {
729: if (!faceGeomFVM) continue;
730: DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
731: field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
732: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
733: }
734: return(0);
735: }
739: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
740: {
741: const PetscInt debug = 0;
742: PetscSection section;
743: PetscQuadrature quad;
744: Vec localX;
745: PetscScalar *funcVal, *interpolant;
746: PetscReal *coords, *v0, *J, *invJ, detJ;
747: PetscReal localDiff = 0.0;
748: const PetscReal *quadPoints, *quadWeights;
749: PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
750: PetscErrorCode ierr;
753: DMGetDimension(dm, &dim);
754: DMGetDefaultSection(dm, §ion);
755: PetscSectionGetNumFields(section, &numFields);
756: DMGetLocalVector(dm, &localX);
757: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
758: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
759: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
760: for (field = 0; field < numFields; ++field) {
761: PetscObject obj;
762: PetscClassId id;
763: PetscInt Nc;
765: DMGetField(dm, field, &obj);
766: PetscObjectGetClassId(obj, &id);
767: if (id == PETSCFE_CLASSID) {
768: PetscFE fe = (PetscFE) obj;
770: PetscFEGetQuadrature(fe, &quad);
771: PetscFEGetNumComponents(fe, &Nc);
772: } else if (id == PETSCFV_CLASSID) {
773: PetscFV fv = (PetscFV) obj;
775: PetscFVGetQuadrature(fv, &quad);
776: PetscFVGetNumComponents(fv, &Nc);
777: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
778: numComponents += Nc;
779: }
780: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
781: PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
782: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
783: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
784: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
785: for (c = cStart; c < cEnd; ++c) {
786: PetscScalar *x = NULL;
787: PetscReal elemDiff = 0.0;
789: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
790: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
791: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
793: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
794: PetscObject obj;
795: PetscClassId id;
796: void * const ctx = ctxs ? ctxs[field] : NULL;
797: PetscInt Nb, Nc, q, fc;
799: DMGetField(dm, field, &obj);
800: PetscObjectGetClassId(obj, &id);
801: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
802: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
803: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
804: if (debug) {
805: char title[1024];
806: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
807: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
808: }
809: for (q = 0; q < numQuadPoints; ++q) {
810: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
811: (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);
812: if (ierr) {
813: PetscErrorCode ierr2;
814: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
815: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
816: ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
817:
818: }
819: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
820: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
821: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
822: for (fc = 0; fc < Nc; ++fc) {
823: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
824: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
825: }
826: }
827: fieldOffset += Nb*Nc;
828: }
829: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
830: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
831: localDiff += elemDiff;
832: }
833: PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
834: DMRestoreLocalVector(dm, &localX);
835: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
836: *diff = PetscSqrtReal(*diff);
837: return(0);
838: }
842: PetscErrorCode DMComputeL2GradientDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
843: {
844: const PetscInt debug = 0;
845: PetscSection section;
846: PetscQuadrature quad;
847: Vec localX;
848: PetscScalar *funcVal, *interpolantVec;
849: PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
850: PetscReal localDiff = 0.0;
851: PetscInt dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
852: PetscErrorCode ierr;
855: DMGetDimension(dm, &dim);
856: DMGetDefaultSection(dm, §ion);
857: PetscSectionGetNumFields(section, &numFields);
858: DMGetLocalVector(dm, &localX);
859: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
860: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
861: for (field = 0; field < numFields; ++field) {
862: PetscFE fe;
863: PetscInt Nc;
865: DMGetField(dm, field, (PetscObject *) &fe);
866: PetscFEGetQuadrature(fe, &quad);
867: PetscFEGetNumComponents(fe, &Nc);
868: numComponents += Nc;
869: }
870: /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
871: PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
872: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
873: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
874: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
875: for (c = cStart; c < cEnd; ++c) {
876: PetscScalar *x = NULL;
877: PetscReal elemDiff = 0.0;
879: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
880: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
881: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
883: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
884: PetscFE fe;
885: void * const ctx = ctxs ? ctxs[field] : NULL;
886: const PetscReal *quadPoints, *quadWeights;
887: PetscReal *basisDer;
888: PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
890: DMGetField(dm, field, (PetscObject *) &fe);
891: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
892: PetscFEGetDimension(fe, &Nb);
893: PetscFEGetNumComponents(fe, &Ncomp);
894: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
895: if (debug) {
896: char title[1024];
897: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
898: DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);
899: }
900: for (q = 0; q < numQuadPoints; ++q) {
901: for (d = 0; d < dim; d++) {
902: coords[d] = v0[d];
903: for (e = 0; e < dim; e++) {
904: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
905: }
906: }
907: (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);
908: if (ierr) {
909: PetscErrorCode ierr2;
910: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
911: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
912: ierr2 = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr2);
913:
914: }
915: for (fc = 0; fc < Ncomp; ++fc) {
916: PetscScalar interpolant = 0.0;
918: for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
919: for (f = 0; f < Nb; ++f) {
920: const PetscInt fidx = f*Ncomp+fc;
922: for (d = 0; d < dim; ++d) {
923: realSpaceDer[d] = 0.0;
924: for (g = 0; g < dim; ++g) {
925: realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
926: }
927: interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
928: }
929: }
930: for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
931: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
932: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
933: }
934: }
935: comp += Ncomp;
936: fieldOffset += Nb*Ncomp;
937: }
938: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
939: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
940: localDiff += elemDiff;
941: }
942: PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
943: DMRestoreLocalVector(dm, &localX);
944: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
945: *diff = PetscSqrtReal(*diff);
946: return(0);
947: }
951: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
952: {
953: const PetscInt debug = 0;
954: PetscSection section;
955: PetscQuadrature quad;
956: Vec localX;
957: PetscScalar *funcVal, *interpolant;
958: PetscReal *coords, *v0, *J, *invJ, detJ;
959: PetscReal *localDiff;
960: const PetscReal *quadPoints, *quadWeights;
961: PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
962: PetscErrorCode ierr;
965: DMGetDimension(dm, &dim);
966: DMGetDefaultSection(dm, §ion);
967: PetscSectionGetNumFields(section, &numFields);
968: DMGetLocalVector(dm, &localX);
969: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
970: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
971: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
972: for (field = 0; field < numFields; ++field) {
973: PetscObject obj;
974: PetscClassId id;
975: PetscInt Nc;
977: DMGetField(dm, field, &obj);
978: PetscObjectGetClassId(obj, &id);
979: if (id == PETSCFE_CLASSID) {
980: PetscFE fe = (PetscFE) obj;
982: PetscFEGetQuadrature(fe, &quad);
983: PetscFEGetNumComponents(fe, &Nc);
984: } else if (id == PETSCFV_CLASSID) {
985: PetscFV fv = (PetscFV) obj;
987: PetscFVGetQuadrature(fv, &quad);
988: PetscFVGetNumComponents(fv, &Nc);
989: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
990: numComponents += Nc;
991: }
992: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
993: PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
994: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
995: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
996: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
997: for (c = cStart; c < cEnd; ++c) {
998: PetscScalar *x = NULL;
1000: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
1001: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1002: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1004: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1005: PetscObject obj;
1006: PetscClassId id;
1007: void * const ctx = ctxs ? ctxs[field] : NULL;
1008: PetscInt Nb, Nc, q, fc;
1010: PetscReal elemDiff = 0.0;
1012: DMGetField(dm, field, &obj);
1013: PetscObjectGetClassId(obj, &id);
1014: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1015: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1016: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1017: if (debug) {
1018: char title[1024];
1019: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1020: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1021: }
1022: for (q = 0; q < numQuadPoints; ++q) {
1023: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1024: (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);
1025: if (ierr) {
1026: PetscErrorCode ierr2;
1027: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1028: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1029: ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
1030:
1031: }
1032: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1033: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1034: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1035: for (fc = 0; fc < Nc; ++fc) {
1036: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
1037: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1038: }
1039: }
1040: fieldOffset += Nb*Nc;
1041: localDiff[field] += elemDiff;
1042: }
1043: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1044: }
1045: DMRestoreLocalVector(dm, &localX);
1046: MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1047: for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1048: PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);
1049: return(0);
1050: }
1054: /*@C
1055: DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec.
1057: Input Parameters:
1058: + dm - The DM
1059: . time - The time
1060: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
1061: . ctxs - Optional array of contexts to pass to each function, or NULL.
1062: - X - The coefficient vector u_h
1064: Output Parameter:
1065: . D - A Vec which holds the difference ||u - u_h||_2 for each cell
1067: Level: developer
1069: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1070: @*/
1071: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
1072: {
1073: PetscSection section;
1074: PetscQuadrature quad;
1075: Vec localX;
1076: PetscScalar *funcVal, *interpolant;
1077: PetscReal *coords, *v0, *J, *invJ, detJ;
1078: const PetscReal *quadPoints, *quadWeights;
1079: PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1080: PetscErrorCode ierr;
1083: VecSet(D, 0.0);
1084: DMGetDimension(dm, &dim);
1085: DMGetDefaultSection(dm, §ion);
1086: PetscSectionGetNumFields(section, &numFields);
1087: DMGetLocalVector(dm, &localX);
1088: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1089: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1090: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1091: for (field = 0; field < numFields; ++field) {
1092: PetscObject obj;
1093: PetscClassId id;
1094: PetscInt Nc;
1096: DMGetField(dm, field, &obj);
1097: PetscObjectGetClassId(obj, &id);
1098: if (id == PETSCFE_CLASSID) {
1099: PetscFE fe = (PetscFE) obj;
1101: PetscFEGetQuadrature(fe, &quad);
1102: PetscFEGetNumComponents(fe, &Nc);
1103: } else if (id == PETSCFV_CLASSID) {
1104: PetscFV fv = (PetscFV) obj;
1106: PetscFVGetQuadrature(fv, &quad);
1107: PetscFVGetNumComponents(fv, &Nc);
1108: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1109: numComponents += Nc;
1110: }
1111: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1112: PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1113: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1114: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1115: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1116: for (c = cStart; c < cEnd; ++c) {
1117: PetscScalar *x = NULL;
1118: PetscScalar elemDiff = 0.0;
1120: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
1121: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1122: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1124: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1125: PetscObject obj;
1126: PetscClassId id;
1127: void * const ctx = ctxs ? ctxs[field] : NULL;
1128: PetscInt Nb, Nc, q, fc;
1130: DMGetField(dm, field, &obj);
1131: PetscObjectGetClassId(obj, &id);
1132: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1133: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1134: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1135: if (funcs[field]) {
1136: for (q = 0; q < numQuadPoints; ++q) {
1137: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1138: (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);
1139: if (ierr) {
1140: PetscErrorCode ierr2;
1141: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1142: ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
1143: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1144:
1145: }
1146: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1147: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1148: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1149: for (fc = 0; fc < Nc; ++fc) {
1150: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1151: }
1152: }
1153: }
1154: fieldOffset += Nb*Nc;
1155: }
1156: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1157: VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1158: }
1159: PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
1160: DMRestoreLocalVector(dm, &localX);
1161: VecSqrtAbs(D);
1162: return(0);
1163: }
1167: /*@
1168: DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user
1170: Input Parameters:
1171: + dm - The mesh
1172: . X - Local input vector
1173: - user - The user context
1175: Output Parameter:
1176: . integral - Local integral for each field
1178: Level: developer
1180: .seealso: DMPlexComputeResidualFEM()
1181: @*/
1182: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1183: {
1184: DM_Plex *mesh = (DM_Plex *) dm->data;
1185: DM dmAux;
1186: Vec localX, A;
1187: PetscDS prob, probAux = NULL;
1188: PetscSection section, sectionAux;
1189: PetscFECellGeom *cgeom;
1190: PetscScalar *u, *a = NULL;
1191: PetscReal *lintegral, *vol;
1192: PetscInt dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c;
1193: PetscInt totDim, totDimAux;
1194: PetscErrorCode ierr;
1197: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1198: DMGetLocalVector(dm, &localX);
1199: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);
1200: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1201: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1202: DMGetDimension(dm, &dim);
1203: DMGetDefaultSection(dm, §ion);
1204: DMGetDS(dm, &prob);
1205: PetscDSGetTotalDimension(prob, &totDim);
1206: PetscSectionGetNumFields(section, &Nf);
1207: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1208: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1209: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1210: numCells = cEnd - cStart;
1211: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1212: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1213: if (dmAux) {
1214: DMGetDefaultSection(dmAux, §ionAux);
1215: DMGetDS(dmAux, &probAux);
1216: PetscDSGetTotalDimension(probAux, &totDimAux);
1217: }
1218: PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);
1219: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1220: for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1221: for (c = cStart; c < cEnd; ++c) {
1222: PetscScalar *x = NULL;
1223: PetscInt i;
1225: DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);
1226: DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);
1227: if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1228: DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1229: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1230: DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1231: if (dmAux) {
1232: DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1233: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1234: DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1235: }
1236: }
1237: for (f = 0; f < Nf; ++f) {
1238: PetscObject obj;
1239: PetscClassId id;
1240: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1242: PetscDSGetDiscretization(prob, f, &obj);
1243: PetscObjectGetClassId(obj, &id);
1244: if (id == PETSCFE_CLASSID) {
1245: PetscFE fe = (PetscFE) obj;
1246: PetscQuadrature q;
1247: PetscInt Nq, Nb;
1249: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1250: PetscFEGetQuadrature(fe, &q);
1251: PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1252: PetscFEGetDimension(fe, &Nb);
1253: blockSize = Nb*Nq;
1254: batchSize = numBlocks * blockSize;
1255: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1256: numChunks = numCells / (numBatches*batchSize);
1257: Ne = numChunks*numBatches*batchSize;
1258: Nr = numCells % (numBatches*batchSize);
1259: offset = numCells - Nr;
1260: PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);
1261: PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1262: } else if (id == PETSCFV_CLASSID) {
1263: /* PetscFV fv = (PetscFV) obj; */
1264: PetscInt foff;
1265: PetscPointFunc obj_func;
1266: PetscScalar lint;
1268: PetscDSGetObjective(prob, f, &obj_func);
1269: PetscDSGetFieldOffset(prob, f, &foff);
1270: if (obj_func) {
1271: for (c = 0; c < numCells; ++c) {
1272: /* TODO: Need full pointwise interpolation and get centroid for x argument */
1273: obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint);
1274: lintegral[f] = PetscRealPart(lint)*vol[c];
1275: }
1276: }
1277: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1278: }
1279: if (dmAux) {PetscFree(a);}
1280: if (mesh->printFEM) {
1281: PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1282: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1283: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1284: }
1285: DMRestoreLocalVector(dm, &localX);
1286: MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1287: PetscFree4(lintegral,u,cgeom,vol);
1288: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1289: return(0);
1290: }
1294: /*@
1295: DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1297: Input Parameters:
1298: + dmf - The fine mesh
1299: . dmc - The coarse mesh
1300: - user - The user context
1302: Output Parameter:
1303: . In - The interpolation matrix
1305: Level: developer
1307: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1308: @*/
1309: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1310: {
1311: DM_Plex *mesh = (DM_Plex *) dmc->data;
1312: const char *name = "Interpolator";
1313: PetscDS prob;
1314: PetscFE *feRef;
1315: PetscFV *fvRef;
1316: PetscSection fsection, fglobalSection;
1317: PetscSection csection, cglobalSection;
1318: PetscScalar *elemMat;
1319: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1320: PetscInt cTotDim, rTotDim = 0;
1321: PetscErrorCode ierr;
1324: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1325: DMGetDimension(dmf, &dim);
1326: DMGetDefaultSection(dmf, &fsection);
1327: DMGetDefaultGlobalSection(dmf, &fglobalSection);
1328: DMGetDefaultSection(dmc, &csection);
1329: DMGetDefaultGlobalSection(dmc, &cglobalSection);
1330: PetscSectionGetNumFields(fsection, &Nf);
1331: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1332: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1333: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1334: DMGetDS(dmf, &prob);
1335: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1336: for (f = 0; f < Nf; ++f) {
1337: PetscObject obj;
1338: PetscClassId id;
1339: PetscInt rNb = 0, Nc = 0;
1341: PetscDSGetDiscretization(prob, f, &obj);
1342: PetscObjectGetClassId(obj, &id);
1343: if (id == PETSCFE_CLASSID) {
1344: PetscFE fe = (PetscFE) obj;
1346: PetscFERefine(fe, &feRef[f]);
1347: PetscFEGetDimension(feRef[f], &rNb);
1348: PetscFEGetNumComponents(fe, &Nc);
1349: } else if (id == PETSCFV_CLASSID) {
1350: PetscFV fv = (PetscFV) obj;
1351: PetscDualSpace Q;
1353: PetscFVRefine(fv, &fvRef[f]);
1354: PetscFVGetDualSpace(fvRef[f], &Q);
1355: PetscDualSpaceGetDimension(Q, &rNb);
1356: PetscFVGetNumComponents(fv, &Nc);
1357: }
1358: rTotDim += rNb*Nc;
1359: }
1360: PetscDSGetTotalDimension(prob, &cTotDim);
1361: PetscMalloc1(rTotDim*cTotDim,&elemMat);
1362: PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1363: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1364: PetscDualSpace Qref;
1365: PetscQuadrature f;
1366: const PetscReal *qpoints, *qweights;
1367: PetscReal *points;
1368: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
1370: /* Compose points from all dual basis functionals */
1371: if (feRef[fieldI]) {
1372: PetscFEGetDualSpace(feRef[fieldI], &Qref);
1373: PetscFEGetNumComponents(feRef[fieldI], &Nc);
1374: } else {
1375: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1376: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1377: }
1378: PetscDualSpaceGetDimension(Qref, &fpdim);
1379: for (i = 0; i < fpdim; ++i) {
1380: PetscDualSpaceGetFunctional(Qref, i, &f);
1381: PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1382: npoints += Np;
1383: }
1384: PetscMalloc1(npoints*dim,&points);
1385: for (i = 0, k = 0; i < fpdim; ++i) {
1386: PetscDualSpaceGetFunctional(Qref, i, &f);
1387: PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1388: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1389: }
1391: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1392: PetscObject obj;
1393: PetscClassId id;
1394: PetscReal *B;
1395: PetscInt NcJ = 0, cpdim = 0, j;
1397: PetscDSGetDiscretization(prob, fieldJ, &obj);
1398: PetscObjectGetClassId(obj, &id);
1399: if (id == PETSCFE_CLASSID) {
1400: PetscFE fe = (PetscFE) obj;
1402: /* Evaluate basis at points */
1403: PetscFEGetNumComponents(fe, &NcJ);
1404: PetscFEGetDimension(fe, &cpdim);
1405: /* For now, fields only interpolate themselves */
1406: if (fieldI == fieldJ) {
1407: 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);
1408: PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1409: for (i = 0, k = 0; i < fpdim; ++i) {
1410: PetscDualSpaceGetFunctional(Qref, i, &f);
1411: PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1412: for (p = 0; p < Np; ++p, ++k) {
1413: for (j = 0; j < cpdim; ++j) {
1414: 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];
1415: }
1416: }
1417: }
1418: PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1419: }
1420: } else if (id == PETSCFV_CLASSID) {
1421: PetscFV fv = (PetscFV) obj;
1423: /* Evaluate constant function at points */
1424: PetscFVGetNumComponents(fv, &NcJ);
1425: cpdim = 1;
1426: /* For now, fields only interpolate themselves */
1427: if (fieldI == fieldJ) {
1428: 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);
1429: for (i = 0, k = 0; i < fpdim; ++i) {
1430: PetscDualSpaceGetFunctional(Qref, i, &f);
1431: PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1432: for (p = 0; p < Np; ++p, ++k) {
1433: for (j = 0; j < cpdim; ++j) {
1434: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1435: }
1436: }
1437: }
1438: }
1439: }
1440: offsetJ += cpdim*NcJ;
1441: }
1442: offsetI += fpdim*Nc;
1443: PetscFree(points);
1444: }
1445: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1446: /* Preallocate matrix */
1447: {
1448: Mat preallocator;
1449: PetscScalar *vals;
1450: PetscInt *cellCIndices, *cellFIndices;
1451: PetscInt locRows, locCols, cell;
1453: MatGetLocalSize(In, &locRows, &locCols);
1454: MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1455: MatSetType(preallocator, MATPREALLOCATOR);
1456: MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1457: MatSetUp(preallocator);
1458: PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1459: for (cell = cStart; cell < cEnd; ++cell) {
1460: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1461: MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1462: }
1463: PetscFree3(vals,cellCIndices,cellFIndices);
1464: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1465: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1466: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1467: MatDestroy(&preallocator);
1468: }
1469: /* Fill matrix */
1470: MatZeroEntries(In);
1471: for (c = cStart; c < cEnd; ++c) {
1472: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1473: }
1474: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1475: PetscFree2(feRef,fvRef);
1476: PetscFree(elemMat);
1477: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1478: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1479: if (mesh->printFEM) {
1480: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1481: MatChop(In, 1.0e-10);
1482: MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1483: }
1484: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1485: return(0);
1486: }
1490: /*@
1491: DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.
1493: Input Parameters:
1494: + dmf - The fine mesh
1495: . dmc - The coarse mesh
1496: - user - The user context
1498: Output Parameter:
1499: . In - The interpolation matrix
1501: Level: developer
1503: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1504: @*/
1505: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1506: {
1507: DM_Plex *mesh = (DM_Plex *) dmf->data;
1508: const char *name = "Interpolator";
1509: PetscDS prob;
1510: PetscSection fsection, csection, globalFSection, globalCSection;
1511: PetscHashJK ht;
1512: PetscLayout rLayout;
1513: PetscInt *dnz, *onz;
1514: PetscInt locRows, rStart, rEnd;
1515: PetscReal *x, *v0, *J, *invJ, detJ;
1516: PetscReal *v0c, *Jc, *invJc, detJc;
1517: PetscScalar *elemMat;
1518: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
1522: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1523: DMGetCoordinateDim(dmc, &dim);
1524: DMGetDS(dmc, &prob);
1525: PetscDSGetRefCoordArrays(prob, &x, NULL);
1526: PetscDSGetNumFields(prob, &Nf);
1527: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1528: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1529: DMGetDefaultSection(dmf, &fsection);
1530: DMGetDefaultGlobalSection(dmf, &globalFSection);
1531: DMGetDefaultSection(dmc, &csection);
1532: DMGetDefaultGlobalSection(dmc, &globalCSection);
1533: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1534: PetscDSGetTotalDimension(prob, &totDim);
1535: PetscMalloc1(totDim*totDim, &elemMat);
1537: MatGetLocalSize(In, &locRows, NULL);
1538: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1539: PetscLayoutSetLocalSize(rLayout, locRows);
1540: PetscLayoutSetBlockSize(rLayout, 1);
1541: PetscLayoutSetUp(rLayout);
1542: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1543: PetscLayoutDestroy(&rLayout);
1544: PetscCalloc2(locRows,&dnz,locRows,&onz);
1545: PetscHashJKCreate(&ht);
1546: for (field = 0; field < Nf; ++field) {
1547: PetscObject obj;
1548: PetscClassId id;
1549: PetscDualSpace Q = NULL;
1550: PetscQuadrature f;
1551: const PetscReal *qpoints;
1552: PetscInt Nc, Np, fpdim, i, d;
1554: PetscDSGetDiscretization(prob, field, &obj);
1555: PetscObjectGetClassId(obj, &id);
1556: if (id == PETSCFE_CLASSID) {
1557: PetscFE fe = (PetscFE) obj;
1559: PetscFEGetDualSpace(fe, &Q);
1560: PetscFEGetNumComponents(fe, &Nc);
1561: } else if (id == PETSCFV_CLASSID) {
1562: PetscFV fv = (PetscFV) obj;
1564: PetscFVGetDualSpace(fv, &Q);
1565: Nc = 1;
1566: }
1567: PetscDualSpaceGetDimension(Q, &fpdim);
1568: /* For each fine grid cell */
1569: for (cell = cStart; cell < cEnd; ++cell) {
1570: PetscInt *findices, *cindices;
1571: PetscInt numFIndices, numCIndices;
1573: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1574: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1575: if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1576: for (i = 0; i < fpdim; ++i) {
1577: Vec pointVec;
1578: PetscScalar *pV;
1579: PetscSF coarseCellSF = NULL;
1580: const PetscSFNode *coarseCells;
1581: PetscInt numCoarseCells, q, r, c;
1583: /* Get points from the dual basis functional quadrature */
1584: PetscDualSpaceGetFunctional(Q, i, &f);
1585: PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1586: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1587: VecSetBlockSize(pointVec, dim);
1588: VecGetArray(pointVec, &pV);
1589: for (q = 0; q < Np; ++q) {
1590: /* Transform point to real space */
1591: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1592: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1593: }
1594: VecRestoreArray(pointVec, &pV);
1595: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1596: DMLocatePoints(dmc, pointVec, &coarseCellSF);
1597: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1598: /* Update preallocation info */
1599: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1600: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1601: for (r = 0; r < Nc; ++r) {
1602: PetscHashJKKey key;
1603: PetscHashJKIter missing, iter;
1605: key.j = findices[i*Nc+r];
1606: if (key.j < 0) continue;
1607: /* Get indices for coarse elements */
1608: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1609: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1610: for (c = 0; c < numCIndices; ++c) {
1611: key.k = cindices[c];
1612: if (key.k < 0) continue;
1613: PetscHashJKPut(ht, key, &missing, &iter);
1614: if (missing) {
1615: PetscHashJKSet(ht, iter, 1);
1616: if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1617: else ++onz[key.j-rStart];
1618: }
1619: }
1620: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1621: }
1622: }
1623: PetscSFDestroy(&coarseCellSF);
1624: VecDestroy(&pointVec);
1625: }
1626: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1627: }
1628: }
1629: PetscHashJKDestroy(&ht);
1630: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1631: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1632: PetscFree2(dnz,onz);
1633: for (field = 0; field < Nf; ++field) {
1634: PetscObject obj;
1635: PetscClassId id;
1636: PetscDualSpace Q = NULL;
1637: PetscQuadrature f;
1638: const PetscReal *qpoints, *qweights;
1639: PetscInt Nc, Np, fpdim, i, d;
1641: PetscDSGetDiscretization(prob, field, &obj);
1642: PetscObjectGetClassId(obj, &id);
1643: if (id == PETSCFE_CLASSID) {
1644: PetscFE fe = (PetscFE) obj;
1646: PetscFEGetDualSpace(fe, &Q);
1647: PetscFEGetNumComponents(fe, &Nc);
1648: } else if (id == PETSCFV_CLASSID) {
1649: PetscFV fv = (PetscFV) obj;
1651: PetscFVGetDualSpace(fv, &Q);
1652: Nc = 1;
1653: }
1654: PetscDualSpaceGetDimension(Q, &fpdim);
1655: /* For each fine grid cell */
1656: for (cell = cStart; cell < cEnd; ++cell) {
1657: PetscInt *findices, *cindices;
1658: PetscInt numFIndices, numCIndices;
1660: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1661: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1662: if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1663: for (i = 0; i < fpdim; ++i) {
1664: Vec pointVec;
1665: PetscScalar *pV;
1666: PetscSF coarseCellSF = NULL;
1667: const PetscSFNode *coarseCells;
1668: PetscInt numCoarseCells, cpdim, q, c, j;
1670: /* Get points from the dual basis functional quadrature */
1671: PetscDualSpaceGetFunctional(Q, i, &f);
1672: PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);
1673: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1674: VecSetBlockSize(pointVec, dim);
1675: VecGetArray(pointVec, &pV);
1676: for (q = 0; q < Np; ++q) {
1677: /* Transform point to real space */
1678: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1679: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1680: }
1681: VecRestoreArray(pointVec, &pV);
1682: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1683: DMLocatePoints(dmc, pointVec, &coarseCellSF);
1684: /* Update preallocation info */
1685: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1686: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1687: VecGetArray(pointVec, &pV);
1688: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1689: PetscReal pVReal[3];
1691: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1692: /* Transform points from real space to coarse reference space */
1693: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
1694: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
1695: CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x);
1697: if (id == PETSCFE_CLASSID) {
1698: PetscFE fe = (PetscFE) obj;
1699: PetscReal *B;
1701: /* Evaluate coarse basis on contained point */
1702: PetscFEGetDimension(fe, &cpdim);
1703: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1704: /* Get elemMat entries by multiplying by weight */
1705: for (j = 0; j < cpdim; ++j) {
1706: for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell];
1707: }
1708: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1709: } else {
1710: cpdim = 1;
1711: for (j = 0; j < cpdim; ++j) {
1712: for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell];
1713: }
1714: }
1715: /* Update interpolator */
1716: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);}
1717: MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);
1718: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1719: }
1720: VecRestoreArray(pointVec, &pV);
1721: PetscSFDestroy(&coarseCellSF);
1722: VecDestroy(&pointVec);
1723: }
1724: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1725: }
1726: }
1727: PetscFree3(v0,J,invJ);
1728: PetscFree3(v0c,Jc,invJc);
1729: PetscFree(elemMat);
1730: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1731: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1732: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1733: return(0);
1734: }
1738: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1739: {
1740: PetscDS prob;
1741: PetscFE *feRef;
1742: PetscFV *fvRef;
1743: Vec fv, cv;
1744: IS fis, cis;
1745: PetscSection fsection, fglobalSection, csection, cglobalSection;
1746: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1747: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
1751: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1752: DMGetDimension(dmf, &dim);
1753: DMGetDefaultSection(dmf, &fsection);
1754: DMGetDefaultGlobalSection(dmf, &fglobalSection);
1755: DMGetDefaultSection(dmc, &csection);
1756: DMGetDefaultGlobalSection(dmc, &cglobalSection);
1757: PetscSectionGetNumFields(fsection, &Nf);
1758: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1759: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1760: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1761: DMGetDS(dmc, &prob);
1762: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1763: for (f = 0; f < Nf; ++f) {
1764: PetscObject obj;
1765: PetscClassId id;
1766: PetscInt fNb = 0, Nc = 0;
1768: PetscDSGetDiscretization(prob, f, &obj);
1769: PetscObjectGetClassId(obj, &id);
1770: if (id == PETSCFE_CLASSID) {
1771: PetscFE fe = (PetscFE) obj;
1773: PetscFERefine(fe, &feRef[f]);
1774: PetscFEGetDimension(feRef[f], &fNb);
1775: PetscFEGetNumComponents(fe, &Nc);
1776: } else if (id == PETSCFV_CLASSID) {
1777: PetscFV fv = (PetscFV) obj;
1778: PetscDualSpace Q;
1780: PetscFVRefine(fv, &fvRef[f]);
1781: PetscFVGetDualSpace(fvRef[f], &Q);
1782: PetscDualSpaceGetDimension(Q, &fNb);
1783: PetscFVGetNumComponents(fv, &Nc);
1784: }
1785: fTotDim += fNb*Nc;
1786: }
1787: PetscDSGetTotalDimension(prob, &cTotDim);
1788: PetscMalloc1(cTotDim,&cmap);
1789: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1790: PetscFE feC;
1791: PetscFV fvC;
1792: PetscDualSpace QF, QC;
1793: PetscInt NcF, NcC, fpdim, cpdim;
1795: if (feRef[field]) {
1796: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1797: PetscFEGetNumComponents(feC, &NcC);
1798: PetscFEGetNumComponents(feRef[field], &NcF);
1799: PetscFEGetDualSpace(feRef[field], &QF);
1800: PetscDualSpaceGetDimension(QF, &fpdim);
1801: PetscFEGetDualSpace(feC, &QC);
1802: PetscDualSpaceGetDimension(QC, &cpdim);
1803: } else {
1804: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1805: PetscFVGetNumComponents(fvC, &NcC);
1806: PetscFVGetNumComponents(fvRef[field], &NcF);
1807: PetscFVGetDualSpace(fvRef[field], &QF);
1808: PetscDualSpaceGetDimension(QF, &fpdim);
1809: PetscFVGetDualSpace(fvC, &QC);
1810: PetscDualSpaceGetDimension(QC, &cpdim);
1811: }
1812: 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);
1813: for (c = 0; c < cpdim; ++c) {
1814: PetscQuadrature cfunc;
1815: const PetscReal *cqpoints;
1816: PetscInt NpC;
1817: PetscBool found = PETSC_FALSE;
1819: PetscDualSpaceGetFunctional(QC, c, &cfunc);
1820: PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1821: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1822: for (f = 0; f < fpdim; ++f) {
1823: PetscQuadrature ffunc;
1824: const PetscReal *fqpoints;
1825: PetscReal sum = 0.0;
1826: PetscInt NpF, comp;
1828: PetscDualSpaceGetFunctional(QF, f, &ffunc);
1829: PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1830: if (NpC != NpF) continue;
1831: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1832: if (sum > 1.0e-9) continue;
1833: for (comp = 0; comp < NcC; ++comp) {
1834: cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1835: }
1836: found = PETSC_TRUE;
1837: break;
1838: }
1839: if (!found) {
1840: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1841: if (fvRef[field]) {
1842: PetscInt comp;
1843: for (comp = 0; comp < NcC; ++comp) {
1844: cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1845: }
1846: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1847: }
1848: }
1849: offsetC += cpdim*NcC;
1850: offsetF += fpdim*NcF;
1851: }
1852: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1853: PetscFree2(feRef,fvRef);
1855: DMGetGlobalVector(dmf, &fv);
1856: DMGetGlobalVector(dmc, &cv);
1857: VecGetOwnershipRange(cv, &startC, &endC);
1858: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1859: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1860: PetscMalloc1(m,&cindices);
1861: PetscMalloc1(m,&findices);
1862: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1863: for (c = cStart; c < cEnd; ++c) {
1864: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1865: for (d = 0; d < cTotDim; ++d) {
1866: if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
1867: 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]]);
1868: cindices[cellCIndices[d]-startC] = cellCIndices[d];
1869: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1870: }
1871: }
1872: PetscFree(cmap);
1873: PetscFree2(cellCIndices,cellFIndices);
1875: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1876: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1877: VecScatterCreate(cv, cis, fv, fis, sc);
1878: ISDestroy(&cis);
1879: ISDestroy(&fis);
1880: DMRestoreGlobalVector(dmf, &fv);
1881: DMRestoreGlobalVector(dmc, &cv);
1882: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1883: return(0);
1884: }