Actual source code: plexfem.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petsc/private/petscfvimpl.h>
7: /*@
8: DMPlexGetScale - Get the scale for the specified fundamental unit
10: Not collective
12: Input Arguments:
13: + dm - the DM
14: - unit - The SI unit
16: Output Argument:
17: . scale - The value used to scale all quantities with this unit
19: Level: advanced
21: .seealso: DMPlexSetScale(), PetscUnit
22: @*/
23: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
24: {
25: DM_Plex *mesh = (DM_Plex*) dm->data;
30: *scale = mesh->scale[unit];
31: return(0);
32: }
34: /*@
35: DMPlexSetScale - Set the scale for the specified fundamental unit
37: Not collective
39: Input Arguments:
40: + dm - the DM
41: . unit - The SI unit
42: - scale - The value used to scale all quantities with this unit
44: Level: advanced
46: .seealso: DMPlexGetScale(), PetscUnit
47: @*/
48: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
49: {
50: DM_Plex *mesh = (DM_Plex*) dm->data;
54: mesh->scale[unit] = scale;
55: return(0);
56: }
58: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
59: {
60: const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
61: PetscInt *ctxInt = (PetscInt *) ctx;
62: PetscInt dim2 = ctxInt[0];
63: PetscInt d = ctxInt[1];
64: PetscInt i, j, k = dim > 2 ? d - dim : d;
67: if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
68: for (i = 0; i < dim; i++) mode[i] = 0.;
69: if (d < dim) {
70: mode[d] = 1.; /* Translation along axis d */
71: } else {
72: for (i = 0; i < dim; i++) {
73: for (j = 0; j < dim; j++) {
74: mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */
75: }
76: }
77: }
78: return(0);
79: }
81: /*@C
82: DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation
84: Collective on DM
86: Input Arguments:
87: . dm - the DM
89: Output Argument:
90: . sp - the null space
92: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
94: Level: advanced
96: .seealso: MatNullSpaceCreate(), PCGAMG
97: @*/
98: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
99: {
100: MPI_Comm comm;
101: Vec mode[6];
102: PetscSection section, globalSection;
103: PetscInt dim, dimEmbed, n, m, d, i, j;
107: PetscObjectGetComm((PetscObject)dm,&comm);
108: DMGetDimension(dm, &dim);
109: DMGetCoordinateDim(dm, &dimEmbed);
110: if (dim == 1) {
111: MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
112: return(0);
113: }
114: DMGetDefaultSection(dm, §ion);
115: DMGetDefaultGlobalSection(dm, &globalSection);
116: PetscSectionGetConstrainedStorageSize(globalSection, &n);
117: m = (dim*(dim+1))/2;
118: VecCreate(comm, &mode[0]);
119: VecSetSizes(mode[0], n, PETSC_DETERMINE);
120: VecSetUp(mode[0]);
121: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
122: for (d = 0; d < m; d++) {
123: PetscInt ctx[2];
124: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
125: void *voidctx = (void *) (&ctx[0]);
127: ctx[0] = dimEmbed;
128: ctx[1] = d;
129: DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
130: }
131: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
132: /* Orthonormalize system */
133: for (i = dim; i < m; ++i) {
134: PetscScalar dots[6];
136: VecMDot(mode[i], i, mode, dots);
137: for (j = 0; j < i; ++j) dots[j] *= -1.0;
138: VecMAXPY(mode[i], i, dots, mode);
139: VecNormalize(mode[i], NULL);
140: }
141: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
142: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
143: return(0);
144: }
146: /*@C
147: DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation
149: Collective on DM
151: Input Arguments:
152: + dm - the DM
153: . nb - The number of bodies
154: . label - The DMLabel marking each domain
155: . nids - The number of ids per body
156: - ids - An array of the label ids in sequence for each domain
158: Output Argument:
159: . sp - the null space
161: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
163: Level: advanced
165: .seealso: MatNullSpaceCreate()
166: @*/
167: PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp)
168: {
169: MPI_Comm comm;
170: PetscSection section, globalSection;
171: Vec *mode;
172: PetscScalar *dots;
173: PetscInt dim, dimEmbed, n, m, b, d, i, j, off;
177: PetscObjectGetComm((PetscObject)dm,&comm);
178: DMGetDimension(dm, &dim);
179: DMGetCoordinateDim(dm, &dimEmbed);
180: DMGetDefaultSection(dm, §ion);
181: DMGetDefaultGlobalSection(dm, &globalSection);
182: PetscSectionGetConstrainedStorageSize(globalSection, &n);
183: m = nb * (dim*(dim+1))/2;
184: PetscMalloc2(m, &mode, m, &dots);
185: VecCreate(comm, &mode[0]);
186: VecSetSizes(mode[0], n, PETSC_DETERMINE);
187: VecSetUp(mode[0]);
188: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
189: for (b = 0, off = 0; b < nb; ++b) {
190: for (d = 0; d < m/nb; ++d) {
191: PetscInt ctx[2];
192: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
193: void *voidctx = (void *) (&ctx[0]);
195: ctx[0] = dimEmbed;
196: ctx[1] = d;
197: DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
198: off += nids[b];
199: }
200: }
201: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
202: /* Orthonormalize system */
203: for (i = 0; i < m; ++i) {
204: VecMDot(mode[i], i, mode, dots);
205: for (j = 0; j < i; ++j) dots[j] *= -1.0;
206: VecMAXPY(mode[i], i, dots, mode);
207: VecNormalize(mode[i], NULL);
208: }
209: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
210: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
211: PetscFree2(mode, dots);
212: return(0);
213: }
215: /*@
216: DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
217: are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
218: evaluating the dual space basis of that point. A basis function is associated with the point in its
219: transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
220: projection height, which is set with this function. By default, the maximum projection height is zero, which means
221: that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior
222: basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
224: Input Parameters:
225: + dm - the DMPlex object
226: - height - the maximum projection height >= 0
228: Level: advanced
230: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
231: @*/
232: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
233: {
234: DM_Plex *plex = (DM_Plex *) dm->data;
238: plex->maxProjectionHeight = height;
239: return(0);
240: }
242: /*@
243: DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
244: DMPlexProjectXXXLocal() functions.
246: Input Parameters:
247: . dm - the DMPlex object
249: Output Parameters:
250: . height - the maximum projection height
252: Level: intermediate
254: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
255: @*/
256: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
257: {
258: DM_Plex *plex = (DM_Plex *) dm->data;
262: *height = plex->maxProjectionHeight;
263: return(0);
264: }
266: /*@C
267: DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector
269: Input Parameters:
270: + dm - The DM, with a PetscDS that matches the problem being constrained
271: . time - The time
272: . field - The field to constrain
273: . Nc - The number of constrained field components, or 0 for all components
274: . comps - An array of constrained component numbers, or NULL for all components
275: . label - The DMLabel defining constrained points
276: . numids - The number of DMLabel ids for constrained points
277: . ids - An array of ids for constrained points
278: . func - A pointwise function giving boundary values
279: - ctx - An optional user context for bcFunc
281: Output Parameter:
282: . locX - A local vector to receives the boundary values
284: Level: developer
286: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
287: @*/
288: PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX)
289: {
290: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
291: void **ctxs;
292: PetscInt numFields;
293: PetscErrorCode ierr;
296: DMGetNumFields(dm, &numFields);
297: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
298: funcs[field] = func;
299: ctxs[field] = ctx;
300: DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
301: PetscFree2(funcs,ctxs);
302: return(0);
303: }
305: /*@C
306: DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector
308: Input Parameters:
309: + dm - The DM, with a PetscDS that matches the problem being constrained
310: . time - The time
311: . locU - A local vector with the input solution values
312: . field - The field to constrain
313: . Nc - The number of constrained field components, or 0 for all components
314: . comps - An array of constrained component numbers, or NULL for all components
315: . label - The DMLabel defining constrained points
316: . numids - The number of DMLabel ids for constrained points
317: . ids - An array of ids for constrained points
318: . func - A pointwise function giving boundary values
319: - ctx - An optional user context for bcFunc
321: Output Parameter:
322: . locX - A local vector to receives the boundary values
324: Level: developer
326: .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary()
327: @*/
328: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
329: void (*func)(PetscInt, PetscInt, PetscInt,
330: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
331: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
332: PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
333: PetscScalar[]),
334: void *ctx, Vec locX)
335: {
336: void (**funcs)(PetscInt, PetscInt, PetscInt,
337: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
338: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
339: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
340: void **ctxs;
341: PetscInt numFields;
342: PetscErrorCode ierr;
345: DMGetNumFields(dm, &numFields);
346: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
347: funcs[field] = func;
348: ctxs[field] = ctx;
349: DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
350: PetscFree2(funcs,ctxs);
351: return(0);
352: }
354: /*@C
355: DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector
357: Input Parameters:
358: + dm - The DM, with a PetscDS that matches the problem being constrained
359: . time - The time
360: . faceGeometry - A vector with the FVM face geometry information
361: . cellGeometry - A vector with the FVM cell geometry information
362: . Grad - A vector with the FVM cell gradient information
363: . field - The field to constrain
364: . Nc - The number of constrained field components, or 0 for all components
365: . comps - An array of constrained component numbers, or NULL for all components
366: . label - The DMLabel defining constrained points
367: . numids - The number of DMLabel ids for constrained points
368: . ids - An array of ids for constrained points
369: . func - A pointwise function giving boundary values
370: - ctx - An optional user context for bcFunc
372: Output Parameter:
373: . locX - A local vector to receives the boundary values
375: Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary()
377: Level: developer
379: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
380: @*/
381: PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
382: PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
383: {
384: PetscDS prob;
385: PetscSF sf;
386: DM dmFace, dmCell, dmGrad;
387: const PetscScalar *facegeom, *cellgeom = NULL, *grad;
388: const PetscInt *leaves;
389: PetscScalar *x, *fx;
390: PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i;
391: PetscErrorCode ierr, ierru = 0;
394: DMGetPointSF(dm, &sf);
395: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
396: nleaves = PetscMax(0, nleaves);
397: DMGetDimension(dm, &dim);
398: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
399: DMGetDS(dm, &prob);
400: VecGetDM(faceGeometry, &dmFace);
401: VecGetArrayRead(faceGeometry, &facegeom);
402: if (cellGeometry) {
403: VecGetDM(cellGeometry, &dmCell);
404: VecGetArrayRead(cellGeometry, &cellgeom);
405: }
406: if (Grad) {
407: PetscFV fv;
409: PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
410: VecGetDM(Grad, &dmGrad);
411: VecGetArrayRead(Grad, &grad);
412: PetscFVGetNumComponents(fv, &pdim);
413: DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
414: }
415: VecGetArray(locX, &x);
416: for (i = 0; i < numids; ++i) {
417: IS faceIS;
418: const PetscInt *faces;
419: PetscInt numFaces, f;
421: DMLabelGetStratumIS(label, ids[i], &faceIS);
422: if (!faceIS) continue; /* No points with that id on this process */
423: ISGetLocalSize(faceIS, &numFaces);
424: ISGetIndices(faceIS, &faces);
425: for (f = 0; f < numFaces; ++f) {
426: const PetscInt face = faces[f], *cells;
427: PetscFVFaceGeom *fg;
429: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
430: PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
431: if (loc >= 0) continue;
432: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
433: DMPlexGetSupport(dm, face, &cells);
434: if (Grad) {
435: PetscFVCellGeom *cg;
436: PetscScalar *cx, *cgrad;
437: PetscScalar *xG;
438: PetscReal dx[3];
439: PetscInt d;
441: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
442: DMPlexPointLocalRead(dm, cells[0], x, &cx);
443: DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
444: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
445: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
446: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
447: ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
448: if (ierru) {
449: ISRestoreIndices(faceIS, &faces);
450: ISDestroy(&faceIS);
451: goto cleanup;
452: }
453: } else {
454: PetscScalar *xI;
455: PetscScalar *xG;
457: DMPlexPointLocalRead(dm, cells[0], x, &xI);
458: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
459: ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
460: if (ierru) {
461: ISRestoreIndices(faceIS, &faces);
462: ISDestroy(&faceIS);
463: goto cleanup;
464: }
465: }
466: }
467: ISRestoreIndices(faceIS, &faces);
468: ISDestroy(&faceIS);
469: }
470: cleanup:
471: VecRestoreArray(locX, &x);
472: if (Grad) {
473: DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
474: VecRestoreArrayRead(Grad, &grad);
475: }
476: if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
477: VecRestoreArrayRead(faceGeometry, &facegeom);
478: CHKERRQ(ierru);
479: return(0);
480: }
482: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
483: {
484: PetscInt numBd, b;
488: PetscDSGetNumBoundary(dm->prob, &numBd);
489: for (b = 0; b < numBd; ++b) {
490: DMBoundaryConditionType type;
491: const char *labelname;
492: DMLabel label;
493: PetscInt field, Nc;
494: const PetscInt *comps;
495: PetscObject obj;
496: PetscClassId id;
497: void (*func)(void);
498: PetscInt numids;
499: const PetscInt *ids;
500: void *ctx;
502: DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
503: if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
504: DMGetLabel(dm, labelname, &label);
505: DMGetField(dm, field, &obj);
506: PetscObjectGetClassId(obj, &id);
507: if (id == PETSCFE_CLASSID) {
508: switch (type) {
509: /* for FEM, there is no insertion to be done for non-essential boundary conditions */
510: case DM_BC_ESSENTIAL:
511: DMPlexLabelAddCells(dm,label);
512: DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
513: DMPlexLabelClearCells(dm,label);
514: break;
515: case DM_BC_ESSENTIAL_FIELD:
516: DMPlexLabelAddCells(dm,label);
517: DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
518: (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
519: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
520: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
521: DMPlexLabelClearCells(dm,label);
522: break;
523: default: break;
524: }
525: } else if (id == PETSCFV_CLASSID) {
526: if (!faceGeomFVM) continue;
527: DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
528: (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
529: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
530: }
531: return(0);
532: }
534: /*@
535: DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector
537: Input Parameters:
538: + dm - The DM
539: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
540: . time - The time
541: . faceGeomFVM - Face geometry data for FV discretizations
542: . cellGeomFVM - Cell geometry data for FV discretizations
543: - gradFVM - Gradient reconstruction data for FV discretizations
545: Output Parameters:
546: . locX - Solution updated with boundary values
548: Level: developer
550: .seealso: DMProjectFunctionLabelLocal()
551: @*/
552: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
553: {
562: PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
563: return(0);
564: }
566: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
567: {
568: Vec localX;
569: PetscErrorCode ierr;
572: DMGetLocalVector(dm, &localX);
573: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
574: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
575: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
576: DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
577: DMRestoreLocalVector(dm, &localX);
578: return(0);
579: }
581: /*@C
582: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
584: Input Parameters:
585: + dm - The DM
586: . time - The time
587: . funcs - The functions to evaluate for each field component
588: . ctxs - Optional array of contexts to pass to each function, or NULL.
589: - localX - The coefficient vector u_h, a local vector
591: Output Parameter:
592: . diff - The diff ||u - u_h||_2
594: Level: developer
596: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
597: @*/
598: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
599: {
600: const PetscInt debug = 0;
601: PetscSection section;
602: PetscQuadrature quad;
603: PetscScalar *funcVal, *interpolant;
604: PetscReal *coords, *detJ, *J;
605: PetscReal localDiff = 0.0;
606: const PetscReal *quadWeights;
607: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset;
608: PetscErrorCode ierr;
611: DMGetDimension(dm, &dim);
612: DMGetCoordinateDim(dm, &coordDim);
613: DMGetDefaultSection(dm, §ion);
614: PetscSectionGetNumFields(section, &numFields);
615: for (field = 0; field < numFields; ++field) {
616: PetscObject obj;
617: PetscClassId id;
618: PetscInt Nc;
620: DMGetField(dm, field, &obj);
621: PetscObjectGetClassId(obj, &id);
622: if (id == PETSCFE_CLASSID) {
623: PetscFE fe = (PetscFE) obj;
625: PetscFEGetQuadrature(fe, &quad);
626: PetscFEGetNumComponents(fe, &Nc);
627: } else if (id == PETSCFV_CLASSID) {
628: PetscFV fv = (PetscFV) obj;
630: PetscFVGetQuadrature(fv, &quad);
631: PetscFVGetNumComponents(fv, &Nc);
632: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
633: numComponents += Nc;
634: }
635: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
636: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
637: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
638: DMPlexGetVTKCellHeight(dm, &cellHeight);
639: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
640: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
641: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
642: for (c = cStart; c < cEnd; ++c) {
643: PetscScalar *x = NULL;
644: PetscReal elemDiff = 0.0;
645: PetscInt qc = 0;
647: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
648: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
650: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
651: PetscObject obj;
652: PetscClassId id;
653: void * const ctx = ctxs ? ctxs[field] : NULL;
654: PetscInt Nb, Nc, q, fc;
656: DMGetField(dm, field, &obj);
657: PetscObjectGetClassId(obj, &id);
658: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
659: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
660: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
661: if (debug) {
662: char title[1024];
663: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
664: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
665: }
666: for (q = 0; q < Nq; ++q) {
667: if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, point %D", detJ[q], c, q);
668: (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
669: if (ierr) {
670: PetscErrorCode ierr2;
671: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
672: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
673: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
674:
675: }
676: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
677: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
678: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
679: for (fc = 0; fc < Nc; ++fc) {
680: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
681: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
682: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
683: }
684: }
685: fieldOffset += Nb;
686: qc += Nc;
687: }
688: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
689: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
690: localDiff += elemDiff;
691: }
692: PetscFree5(funcVal,interpolant,coords,detJ,J);
693: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
694: *diff = PetscSqrtReal(*diff);
695: return(0);
696: }
698: 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)
699: {
700: const PetscInt debug = 0;
701: PetscSection section;
702: PetscQuadrature quad;
703: Vec localX;
704: PetscScalar *funcVal, *interpolant;
705: const PetscReal *quadPoints, *quadWeights;
706: PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ;
707: PetscReal localDiff = 0.0;
708: PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
709: PetscErrorCode ierr;
712: DMGetDimension(dm, &dim);
713: DMGetCoordinateDim(dm, &coordDim);
714: DMGetDefaultSection(dm, §ion);
715: PetscSectionGetNumFields(section, &numFields);
716: DMGetLocalVector(dm, &localX);
717: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
718: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
719: for (field = 0; field < numFields; ++field) {
720: PetscFE fe;
721: PetscInt Nc;
723: DMGetField(dm, field, (PetscObject *) &fe);
724: PetscFEGetQuadrature(fe, &quad);
725: PetscFEGetNumComponents(fe, &Nc);
726: numComponents += Nc;
727: }
728: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
729: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
730: /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
731: PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
732: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
733: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
734: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
735: for (c = cStart; c < cEnd; ++c) {
736: PetscScalar *x = NULL;
737: PetscReal elemDiff = 0.0;
738: PetscInt qc = 0;
740: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);
741: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
743: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
744: PetscFE fe;
745: void * const ctx = ctxs ? ctxs[field] : NULL;
746: PetscReal *basisDer;
747: PetscInt Nb, Nc, q, fc;
749: DMGetField(dm, field, (PetscObject *) &fe);
750: PetscFEGetDimension(fe, &Nb);
751: PetscFEGetNumComponents(fe, &Nc);
752: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
753: if (debug) {
754: char title[1024];
755: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
756: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
757: }
758: for (q = 0; q < Nq; ++q) {
759: if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", detJ[q], c, q);
760: (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
761: if (ierr) {
762: PetscErrorCode ierr2;
763: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
764: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
765: ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
766:
767: }
768: PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
769: for (fc = 0; fc < Nc; ++fc) {
770: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
771: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
772: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
773: }
774: }
775: fieldOffset += Nb;
776: qc += Nc;
777: }
778: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
779: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
780: localDiff += elemDiff;
781: }
782: PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
783: DMRestoreLocalVector(dm, &localX);
784: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
785: *diff = PetscSqrtReal(*diff);
786: return(0);
787: }
789: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
790: {
791: const PetscInt debug = 0;
792: PetscSection section;
793: PetscQuadrature quad;
794: Vec localX;
795: PetscScalar *funcVal, *interpolant;
796: PetscReal *coords, *detJ, *J;
797: PetscReal *localDiff;
798: const PetscReal *quadPoints, *quadWeights;
799: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
800: PetscErrorCode ierr;
803: DMGetDimension(dm, &dim);
804: DMGetCoordinateDim(dm, &coordDim);
805: DMGetDefaultSection(dm, §ion);
806: PetscSectionGetNumFields(section, &numFields);
807: DMGetLocalVector(dm, &localX);
808: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
809: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
810: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
811: for (field = 0; field < numFields; ++field) {
812: PetscObject obj;
813: PetscClassId id;
814: PetscInt Nc;
816: DMGetField(dm, field, &obj);
817: PetscObjectGetClassId(obj, &id);
818: if (id == PETSCFE_CLASSID) {
819: PetscFE fe = (PetscFE) obj;
821: PetscFEGetQuadrature(fe, &quad);
822: PetscFEGetNumComponents(fe, &Nc);
823: } else if (id == PETSCFV_CLASSID) {
824: PetscFV fv = (PetscFV) obj;
826: PetscFVGetQuadrature(fv, &quad);
827: PetscFVGetNumComponents(fv, &Nc);
828: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
829: numComponents += Nc;
830: }
831: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
832: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
833: PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
834: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
835: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
836: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
837: for (c = cStart; c < cEnd; ++c) {
838: PetscScalar *x = NULL;
839: PetscInt qc = 0;
841: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
842: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
844: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
845: PetscObject obj;
846: PetscClassId id;
847: void * const ctx = ctxs ? ctxs[field] : NULL;
848: PetscInt Nb, Nc, q, fc;
850: PetscReal elemDiff = 0.0;
852: DMGetField(dm, field, &obj);
853: PetscObjectGetClassId(obj, &id);
854: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
855: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
856: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
857: if (debug) {
858: char title[1024];
859: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
860: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
861: }
862: for (q = 0; q < Nq; ++q) {
863: if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature point %D", detJ, c, q);
864: (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
865: if (ierr) {
866: PetscErrorCode ierr2;
867: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
868: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
869: ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
870:
871: }
872: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
873: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
874: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
875: for (fc = 0; fc < Nc; ++fc) {
876: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
877: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d point %g %g %g diff %g\n", c, field, coordDim > 0 ? coords[0] : 0., coordDim > 1 ? coords[1] : 0., coordDim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
878: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
879: }
880: }
881: fieldOffset += Nb;
882: qc += Nc;
883: localDiff[field] += elemDiff;
884: }
885: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
886: }
887: DMRestoreLocalVector(dm, &localX);
888: MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
889: for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
890: PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
891: return(0);
892: }
894: /*@C
895: 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.
897: Input Parameters:
898: + dm - The DM
899: . time - The time
900: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
901: . ctxs - Optional array of contexts to pass to each function, or NULL.
902: - X - The coefficient vector u_h
904: Output Parameter:
905: . D - A Vec which holds the difference ||u - u_h||_2 for each cell
907: Level: developer
909: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
910: @*/
911: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
912: {
913: PetscSection section;
914: PetscQuadrature quad;
915: Vec localX;
916: PetscScalar *funcVal, *interpolant;
917: PetscReal *coords, *detJ, *J;
918: const PetscReal *quadPoints, *quadWeights;
919: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
920: PetscErrorCode ierr;
923: VecSet(D, 0.0);
924: DMGetDimension(dm, &dim);
925: DMGetCoordinateDim(dm, &coordDim);
926: DMGetDefaultSection(dm, §ion);
927: PetscSectionGetNumFields(section, &numFields);
928: DMGetLocalVector(dm, &localX);
929: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
930: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
931: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
932: for (field = 0; field < numFields; ++field) {
933: PetscObject obj;
934: PetscClassId id;
935: PetscInt Nc;
937: DMGetField(dm, field, &obj);
938: PetscObjectGetClassId(obj, &id);
939: if (id == PETSCFE_CLASSID) {
940: PetscFE fe = (PetscFE) obj;
942: PetscFEGetQuadrature(fe, &quad);
943: PetscFEGetNumComponents(fe, &Nc);
944: } else if (id == PETSCFV_CLASSID) {
945: PetscFV fv = (PetscFV) obj;
947: PetscFVGetQuadrature(fv, &quad);
948: PetscFVGetNumComponents(fv, &Nc);
949: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
950: numComponents += Nc;
951: }
952: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
953: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
954: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
955: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
956: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
957: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
958: for (c = cStart; c < cEnd; ++c) {
959: PetscScalar *x = NULL;
960: PetscScalar elemDiff = 0.0;
961: PetscInt qc = 0;
963: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
964: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
966: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
967: PetscObject obj;
968: PetscClassId id;
969: void * const ctx = ctxs ? ctxs[field] : NULL;
970: PetscInt Nb, Nc, q, fc;
972: DMGetField(dm, field, &obj);
973: PetscObjectGetClassId(obj, &id);
974: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
975: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
976: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
977: if (funcs[field]) {
978: for (q = 0; q < Nq; ++q) {
979: if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], c, q);
980: (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx);
981: if (ierr) {
982: PetscErrorCode ierr2;
983: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
984: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
985: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
986:
987: }
988: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
989: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
990: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
991: for (fc = 0; fc < Nc; ++fc) {
992: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
993: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
994: }
995: }
996: }
997: fieldOffset += Nb;
998: qc += Nc;
999: }
1000: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1001: VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1002: }
1003: PetscFree5(funcVal,interpolant,coords,detJ,J);
1004: DMRestoreLocalVector(dm, &localX);
1005: VecSqrtAbs(D);
1006: return(0);
1007: }
1009: /*@C
1010: DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec.
1012: Input Parameters:
1013: + dm - The DM
1014: - LocX - The coefficient vector u_h
1016: Output Parameter:
1017: . locC - A Vec which holds the Clement interpolant of the gradient
1019: Notes: Add citation to (Clement, 1975) and definition of the interpolant
1020: \nabla u_h(v_i) = \sum_{T_i \in support(v_i)} |T_i| \nabla u_h(T_i) / \sum_{T_i \in support(v_i)} |T_i| where |T_i| is the cell volume
1022: Level: developer
1024: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1025: @*/
1026: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1027: {
1028: DM_Plex *mesh = (DM_Plex *) dm->data;
1029: PetscInt debug = mesh->printFEM;
1030: DM dmC;
1031: PetscSection section;
1032: PetscQuadrature quad;
1033: PetscScalar *interpolant, *gradsum;
1034: PetscReal *coords, *detJ, *J, *invJ;
1035: const PetscReal *quadPoints, *quadWeights;
1036: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset;
1037: PetscErrorCode ierr;
1040: VecGetDM(locC, &dmC);
1041: VecSet(locC, 0.0);
1042: DMGetDimension(dm, &dim);
1043: DMGetCoordinateDim(dm, &coordDim);
1044: DMGetDefaultSection(dm, §ion);
1045: PetscSectionGetNumFields(section, &numFields);
1046: for (field = 0; field < numFields; ++field) {
1047: PetscObject obj;
1048: PetscClassId id;
1049: PetscInt Nc;
1051: DMGetField(dm, field, &obj);
1052: PetscObjectGetClassId(obj, &id);
1053: if (id == PETSCFE_CLASSID) {
1054: PetscFE fe = (PetscFE) obj;
1056: PetscFEGetQuadrature(fe, &quad);
1057: PetscFEGetNumComponents(fe, &Nc);
1058: } else if (id == PETSCFV_CLASSID) {
1059: PetscFV fv = (PetscFV) obj;
1061: PetscFVGetQuadrature(fv, &quad);
1062: PetscFVGetNumComponents(fv, &Nc);
1063: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1064: numComponents += Nc;
1065: }
1066: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1067: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1068: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);
1069: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1070: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1071: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1072: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1073: for (v = vStart; v < vEnd; ++v) {
1074: PetscScalar volsum = 0.0;
1075: PetscInt *star = NULL;
1076: PetscInt starSize, st, d, fc;
1078: PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));
1079: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1080: for (st = 0; st < starSize*2; st += 2) {
1081: const PetscInt cell = star[st];
1082: PetscScalar *grad = &gradsum[coordDim*numComponents];
1083: PetscScalar *x = NULL;
1084: PetscReal vol = 0.0;
1086: if ((cell < cStart) || (cell >= cEnd)) continue;
1087: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);
1088: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1089: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1090: PetscObject obj;
1091: PetscClassId id;
1092: PetscInt Nb, Nc, q, qc = 0;
1094: PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));
1095: DMGetField(dm, field, &obj);
1096: PetscObjectGetClassId(obj, &id);
1097: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1098: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1099: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1100: for (q = 0; q < Nq; ++q) {
1101: if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], cell, q);
1102: if (ierr) {
1103: PetscErrorCode ierr2;
1104: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1105: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1106: ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2);
1107:
1108: }
1109: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);}
1110: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1111: for (fc = 0; fc < Nc; ++fc) {
1112: const PetscReal wt = quadWeights[q*qNc+qc+fc];
1114: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q];
1115: }
1116: vol += quadWeights[q*qNc]*detJ[q];
1117: }
1118: fieldOffset += Nb;
1119: qc += Nc;
1120: }
1121: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1122: for (fc = 0; fc < numComponents; ++fc) {
1123: for (d = 0; d < coordDim; ++d) {
1124: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1125: }
1126: }
1127: volsum += vol;
1128: if (debug) {
1129: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1130: for (fc = 0; fc < numComponents; ++fc) {
1131: for (d = 0; d < coordDim; ++d) {
1132: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1133: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1134: }
1135: }
1136: PetscPrintf(PETSC_COMM_SELF, "]\n");
1137: }
1138: }
1139: for (fc = 0; fc < numComponents; ++fc) {
1140: for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1141: }
1142: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1143: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1144: }
1145: PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);
1146: return(0);
1147: }
1149: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1150: {
1151: DM dmAux = NULL;
1152: PetscDS prob, probAux;
1153: PetscSection section, sectionAux;
1154: Vec locX, locA;
1155: PetscInt dim, numCells = cEnd - cStart, cell, c, f;
1156: PetscBool useFEM = PETSC_FALSE, useFVM = PETSC_FALSE;
1157: /* DS */
1158: PetscInt Nf, totDim, *uOff, *uOff_x, numConstants;
1159: PetscInt NfAux, totDimAux, *aOff;
1160: PetscScalar *u, *a;
1161: const PetscScalar *constants;
1162: /* Geometry */
1163: PetscFECellGeom *cgeomFEM;
1164: DM dmGrad;
1165: Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1166: PetscFVCellGeom *cgeomFVM;
1167: const PetscScalar *lgrad;
1168: PetscErrorCode ierr;
1171: DMGetDS(dm, &prob);
1172: DMGetDimension(dm, &dim);
1173: DMGetDefaultSection(dm, §ion);
1174: PetscSectionGetNumFields(section, &Nf);
1175: /* Determine which discretizations we have */
1176: for (f = 0; f < Nf; ++f) {
1177: PetscObject obj;
1178: PetscClassId id;
1180: PetscDSGetDiscretization(prob, f, &obj);
1181: PetscObjectGetClassId(obj, &id);
1182: if (id == PETSCFE_CLASSID) useFEM = PETSC_TRUE;
1183: if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1184: }
1185: /* Get local solution with boundary values */
1186: DMGetLocalVector(dm, &locX);
1187: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1188: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1189: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1190: /* Read DS information */
1191: PetscDSGetTotalDimension(prob, &totDim);
1192: PetscDSGetComponentOffsets(prob, &uOff);
1193: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1194: PetscDSGetConstants(prob, &numConstants, &constants);
1195: /* Read Auxiliary DS information */
1196: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1197: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1198: if (dmAux) {
1199: DMGetDS(dmAux, &probAux);
1200: PetscDSGetNumFields(probAux, &NfAux);
1201: DMGetDefaultSection(dmAux, §ionAux);
1202: PetscDSGetTotalDimension(probAux, &totDimAux);
1203: PetscDSGetComponentOffsets(probAux, &aOff);
1204: }
1205: /* Allocate data arrays */
1206: PetscCalloc1(numCells*totDim, &u);
1207: if (useFEM) {PetscMalloc1(numCells, &cgeomFEM);}
1208: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1209: /* Read out geometry */
1210: if (useFEM) {
1211: for (cell = cStart; cell < cEnd; ++cell) {
1212: const PetscInt c = cell - cStart;
1214: DMPlexComputeCellGeometryFEM(dm, cell, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);
1215: if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c);
1216: }
1217: }
1218: if (useFVM) {
1219: PetscFV fv = NULL;
1220: Vec grad;
1221: PetscInt fStart, fEnd;
1222: PetscBool compGrad;
1224: for (f = 0; f < Nf; ++f) {
1225: PetscObject obj;
1226: PetscClassId id;
1228: PetscDSGetDiscretization(prob, f, &obj);
1229: PetscObjectGetClassId(obj, &id);
1230: if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1231: }
1232: PetscFVGetComputeGradients(fv, &compGrad);
1233: PetscFVSetComputeGradients(fv, PETSC_TRUE);
1234: DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1235: DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1236: PetscFVSetComputeGradients(fv, compGrad);
1237: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1238: /* Reconstruct and limit cell gradients */
1239: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1240: DMGetGlobalVector(dmGrad, &grad);
1241: DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1242: /* Communicate gradient values */
1243: DMGetLocalVector(dmGrad, &locGrad);
1244: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1245: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1246: DMRestoreGlobalVector(dmGrad, &grad);
1247: /* Handle non-essential (e.g. outflow) boundary values */
1248: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1249: VecGetArrayRead(locGrad, &lgrad);
1250: }
1251: /* Read out data from inputs */
1252: for (c = cStart; c < cEnd; ++c) {
1253: PetscScalar *x = NULL;
1254: PetscInt i;
1256: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1257: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1258: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1259: if (dmAux) {
1260: DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1261: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1262: DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1263: }
1264: }
1265: /* Do integration for each field */
1266: for (f = 0; f < Nf; ++f) {
1267: PetscObject obj;
1268: PetscClassId id;
1269: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1271: PetscDSGetDiscretization(prob, f, &obj);
1272: PetscObjectGetClassId(obj, &id);
1273: if (id == PETSCFE_CLASSID) {
1274: PetscFE fe = (PetscFE) obj;
1275: PetscQuadrature q;
1276: PetscInt Nq, Nb;
1278: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1279: PetscFEGetQuadrature(fe, &q);
1280: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1281: PetscFEGetDimension(fe, &Nb);
1282: blockSize = Nb*Nq;
1283: batchSize = numBlocks * blockSize;
1284: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1285: numChunks = numCells / (numBatches*batchSize);
1286: Ne = numChunks*numBatches*batchSize;
1287: Nr = numCells % (numBatches*batchSize);
1288: offset = numCells - Nr;
1289: PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, cintegral);
1290: PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1291: } else if (id == PETSCFV_CLASSID) {
1292: PetscInt foff;
1293: PetscPointFunc obj_func;
1294: PetscScalar lint;
1296: PetscDSGetObjective(prob, f, &obj_func);
1297: PetscDSGetFieldOffset(prob, f, &foff);
1298: if (obj_func) {
1299: for (c = 0; c < numCells; ++c) {
1300: PetscScalar *u_x;
1302: DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1303: obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, &a[totDimAux*c], NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint);
1304: cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1305: }
1306: }
1307: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1308: }
1309: /* Cleanup data arrays */
1310: if (useFVM) {
1311: VecRestoreArrayRead(locGrad, &lgrad);
1312: VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1313: DMRestoreLocalVector(dmGrad, &locGrad);
1314: VecDestroy(&faceGeometryFVM);
1315: VecDestroy(&cellGeometryFVM);
1316: DMDestroy(&dmGrad);
1317: }
1318: if (useFEM) {PetscFree(cgeomFEM);}
1319: if (dmAux) {PetscFree(a);}
1320: PetscFree(u);
1321: /* Cleanup */
1322: DMRestoreLocalVector(dm, &locX);
1323: return(0);
1324: }
1326: /*@
1327: DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user
1329: Input Parameters:
1330: + dm - The mesh
1331: . X - Global input vector
1332: - user - The user context
1334: Output Parameter:
1335: . integral - Integral for each field
1337: Level: developer
1339: .seealso: DMPlexComputeResidualFEM()
1340: @*/
1341: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1342: {
1343: DM_Plex *mesh = (DM_Plex *) dm->data;
1344: PetscScalar *cintegral, *lintegral;
1345: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1352: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1353: DMGetNumFields(dm, &Nf);
1354: DMPlexGetVTKCellHeight(dm, &cellHeight);
1355: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1356: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1357: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1358: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1359: PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1360: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1361: /* Sum up values */
1362: for (cell = cStart; cell < cEnd; ++cell) {
1363: const PetscInt c = cell - cStart;
1365: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1366: for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1367: }
1368: MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1369: if (mesh->printFEM) {
1370: PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1371: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1372: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1373: }
1374: PetscFree2(lintegral, cintegral);
1375: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1376: return(0);
1377: }
1379: /*@
1380: DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user
1382: Input Parameters:
1383: + dm - The mesh
1384: . X - Global input vector
1385: - user - The user context
1387: Output Parameter:
1388: . integral - Cellwise integrals for each field
1390: Level: developer
1392: .seealso: DMPlexComputeResidualFEM()
1393: @*/
1394: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
1395: {
1396: DM_Plex *mesh = (DM_Plex *) dm->data;
1397: DM dmF;
1398: PetscSection sectionF;
1399: PetscScalar *cintegral, *af;
1400: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1407: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1408: DMGetNumFields(dm, &Nf);
1409: DMPlexGetVTKCellHeight(dm, &cellHeight);
1410: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1411: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1412: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1413: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1414: PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
1415: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1416: /* Put values in F*/
1417: VecGetDM(F, &dmF);
1418: DMGetDefaultSection(dmF, §ionF);
1419: VecGetArray(F, &af);
1420: for (cell = cStart; cell < cEnd; ++cell) {
1421: const PetscInt c = cell - cStart;
1422: PetscInt dof, off;
1424: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1425: PetscSectionGetDof(sectionF, cell, &dof);
1426: PetscSectionGetOffset(sectionF, cell, &off);
1427: if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
1428: for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
1429: }
1430: VecRestoreArray(F, &af);
1431: PetscFree(cintegral);
1432: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1433: return(0);
1434: }
1436: /*@
1437: DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1439: Input Parameters:
1440: + dmf - The fine mesh
1441: . dmc - The coarse mesh
1442: - user - The user context
1444: Output Parameter:
1445: . In - The interpolation matrix
1447: Level: developer
1449: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1450: @*/
1451: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1452: {
1453: DM_Plex *mesh = (DM_Plex *) dmc->data;
1454: const char *name = "Interpolator";
1455: PetscDS prob;
1456: PetscFE *feRef;
1457: PetscFV *fvRef;
1458: PetscSection fsection, fglobalSection;
1459: PetscSection csection, cglobalSection;
1460: PetscScalar *elemMat;
1461: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1462: PetscInt cTotDim, rTotDim = 0;
1463: PetscErrorCode ierr;
1466: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1467: DMGetDimension(dmf, &dim);
1468: DMGetDefaultSection(dmf, &fsection);
1469: DMGetDefaultGlobalSection(dmf, &fglobalSection);
1470: DMGetDefaultSection(dmc, &csection);
1471: DMGetDefaultGlobalSection(dmc, &cglobalSection);
1472: PetscSectionGetNumFields(fsection, &Nf);
1473: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1474: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1475: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1476: DMGetDS(dmf, &prob);
1477: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1478: for (f = 0; f < Nf; ++f) {
1479: PetscObject obj;
1480: PetscClassId id;
1481: PetscInt rNb = 0, Nc = 0;
1483: PetscDSGetDiscretization(prob, f, &obj);
1484: PetscObjectGetClassId(obj, &id);
1485: if (id == PETSCFE_CLASSID) {
1486: PetscFE fe = (PetscFE) obj;
1488: PetscFERefine(fe, &feRef[f]);
1489: PetscFEGetDimension(feRef[f], &rNb);
1490: PetscFEGetNumComponents(fe, &Nc);
1491: } else if (id == PETSCFV_CLASSID) {
1492: PetscFV fv = (PetscFV) obj;
1493: PetscDualSpace Q;
1495: PetscFVRefine(fv, &fvRef[f]);
1496: PetscFVGetDualSpace(fvRef[f], &Q);
1497: PetscDualSpaceGetDimension(Q, &rNb);
1498: PetscFVGetNumComponents(fv, &Nc);
1499: }
1500: rTotDim += rNb;
1501: }
1502: PetscDSGetTotalDimension(prob, &cTotDim);
1503: PetscMalloc1(rTotDim*cTotDim,&elemMat);
1504: PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1505: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1506: PetscDualSpace Qref;
1507: PetscQuadrature f;
1508: const PetscReal *qpoints, *qweights;
1509: PetscReal *points;
1510: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
1512: /* Compose points from all dual basis functionals */
1513: if (feRef[fieldI]) {
1514: PetscFEGetDualSpace(feRef[fieldI], &Qref);
1515: PetscFEGetNumComponents(feRef[fieldI], &Nc);
1516: } else {
1517: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1518: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1519: }
1520: PetscDualSpaceGetDimension(Qref, &fpdim);
1521: for (i = 0; i < fpdim; ++i) {
1522: PetscDualSpaceGetFunctional(Qref, i, &f);
1523: PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1524: npoints += Np;
1525: }
1526: PetscMalloc1(npoints*dim,&points);
1527: for (i = 0, k = 0; i < fpdim; ++i) {
1528: PetscDualSpaceGetFunctional(Qref, i, &f);
1529: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1530: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1531: }
1533: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1534: PetscObject obj;
1535: PetscClassId id;
1536: PetscReal *B;
1537: PetscInt NcJ = 0, cpdim = 0, j, qNc;
1539: PetscDSGetDiscretization(prob, fieldJ, &obj);
1540: PetscObjectGetClassId(obj, &id);
1541: if (id == PETSCFE_CLASSID) {
1542: PetscFE fe = (PetscFE) obj;
1544: /* Evaluate basis at points */
1545: PetscFEGetNumComponents(fe, &NcJ);
1546: PetscFEGetDimension(fe, &cpdim);
1547: /* For now, fields only interpolate themselves */
1548: if (fieldI == fieldJ) {
1549: 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);
1550: PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1551: for (i = 0, k = 0; i < fpdim; ++i) {
1552: PetscDualSpaceGetFunctional(Qref, i, &f);
1553: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1554: if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
1555: for (p = 0; p < Np; ++p, ++k) {
1556: for (j = 0; j < cpdim; ++j) {
1557: /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */
1558: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1559: }
1560: }
1561: }
1562: PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1563: }
1564: } else if (id == PETSCFV_CLASSID) {
1565: PetscFV fv = (PetscFV) obj;
1567: /* Evaluate constant function at points */
1568: PetscFVGetNumComponents(fv, &NcJ);
1569: cpdim = 1;
1570: /* For now, fields only interpolate themselves */
1571: if (fieldI == fieldJ) {
1572: 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);
1573: for (i = 0, k = 0; i < fpdim; ++i) {
1574: PetscDualSpaceGetFunctional(Qref, i, &f);
1575: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1576: if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
1577: for (p = 0; p < Np; ++p, ++k) {
1578: for (j = 0; j < cpdim; ++j) {
1579: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c];
1580: }
1581: }
1582: }
1583: }
1584: }
1585: offsetJ += cpdim*NcJ;
1586: }
1587: offsetI += fpdim*Nc;
1588: PetscFree(points);
1589: }
1590: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1591: /* Preallocate matrix */
1592: {
1593: Mat preallocator;
1594: PetscScalar *vals;
1595: PetscInt *cellCIndices, *cellFIndices;
1596: PetscInt locRows, locCols, cell;
1598: MatGetLocalSize(In, &locRows, &locCols);
1599: MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1600: MatSetType(preallocator, MATPREALLOCATOR);
1601: MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1602: MatSetUp(preallocator);
1603: PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1604: for (cell = cStart; cell < cEnd; ++cell) {
1605: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1606: MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1607: }
1608: PetscFree3(vals,cellCIndices,cellFIndices);
1609: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1610: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1611: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1612: MatDestroy(&preallocator);
1613: }
1614: /* Fill matrix */
1615: MatZeroEntries(In);
1616: for (c = cStart; c < cEnd; ++c) {
1617: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1618: }
1619: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1620: PetscFree2(feRef,fvRef);
1621: PetscFree(elemMat);
1622: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1623: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1624: if (mesh->printFEM) {
1625: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1626: MatChop(In, 1.0e-10);
1627: MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1628: }
1629: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1630: return(0);
1631: }
1633: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
1634: {
1635: SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
1636: }
1638: /*@
1639: DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.
1641: Input Parameters:
1642: + dmf - The fine mesh
1643: . dmc - The coarse mesh
1644: - user - The user context
1646: Output Parameter:
1647: . In - The interpolation matrix
1649: Level: developer
1651: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1652: @*/
1653: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1654: {
1655: DM_Plex *mesh = (DM_Plex *) dmf->data;
1656: const char *name = "Interpolator";
1657: PetscDS prob;
1658: PetscSection fsection, csection, globalFSection, globalCSection;
1659: PetscHashJK ht;
1660: PetscLayout rLayout;
1661: PetscInt *dnz, *onz;
1662: PetscInt locRows, rStart, rEnd;
1663: PetscReal *x, *v0, *J, *invJ, detJ;
1664: PetscReal *v0c, *Jc, *invJc, detJc;
1665: PetscScalar *elemMat;
1666: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
1670: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1671: DMGetCoordinateDim(dmc, &dim);
1672: DMGetDS(dmc, &prob);
1673: PetscDSGetRefCoordArrays(prob, &x, NULL);
1674: PetscDSGetNumFields(prob, &Nf);
1675: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1676: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1677: DMGetDefaultSection(dmf, &fsection);
1678: DMGetDefaultGlobalSection(dmf, &globalFSection);
1679: DMGetDefaultSection(dmc, &csection);
1680: DMGetDefaultGlobalSection(dmc, &globalCSection);
1681: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1682: PetscDSGetTotalDimension(prob, &totDim);
1683: PetscMalloc1(totDim, &elemMat);
1685: MatGetLocalSize(In, &locRows, NULL);
1686: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1687: PetscLayoutSetLocalSize(rLayout, locRows);
1688: PetscLayoutSetBlockSize(rLayout, 1);
1689: PetscLayoutSetUp(rLayout);
1690: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1691: PetscLayoutDestroy(&rLayout);
1692: PetscCalloc2(locRows,&dnz,locRows,&onz);
1693: PetscHashJKCreate(&ht);
1694: for (field = 0; field < Nf; ++field) {
1695: PetscObject obj;
1696: PetscClassId id;
1697: PetscDualSpace Q = NULL;
1698: PetscQuadrature f;
1699: const PetscReal *qpoints;
1700: PetscInt Nc, Np, fpdim, i, d;
1702: PetscDSGetDiscretization(prob, field, &obj);
1703: PetscObjectGetClassId(obj, &id);
1704: if (id == PETSCFE_CLASSID) {
1705: PetscFE fe = (PetscFE) obj;
1707: PetscFEGetDualSpace(fe, &Q);
1708: PetscFEGetNumComponents(fe, &Nc);
1709: } else if (id == PETSCFV_CLASSID) {
1710: PetscFV fv = (PetscFV) obj;
1712: PetscFVGetDualSpace(fv, &Q);
1713: Nc = 1;
1714: }
1715: PetscDualSpaceGetDimension(Q, &fpdim);
1716: /* For each fine grid cell */
1717: for (cell = cStart; cell < cEnd; ++cell) {
1718: PetscInt *findices, *cindices;
1719: PetscInt numFIndices, numCIndices;
1721: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1722: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1723: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1724: for (i = 0; i < fpdim; ++i) {
1725: Vec pointVec;
1726: PetscScalar *pV;
1727: PetscSF coarseCellSF = NULL;
1728: const PetscSFNode *coarseCells;
1729: PetscInt numCoarseCells, q, c;
1731: /* Get points from the dual basis functional quadrature */
1732: PetscDualSpaceGetFunctional(Q, i, &f);
1733: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1734: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1735: VecSetBlockSize(pointVec, dim);
1736: VecGetArray(pointVec, &pV);
1737: for (q = 0; q < Np; ++q) {
1738: /* Transform point to real space */
1739: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1740: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1741: }
1742: VecRestoreArray(pointVec, &pV);
1743: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1744: /* OPT: Pack all quad points from fine cell */
1745: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1746: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1747: /* Update preallocation info */
1748: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1749: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1750: {
1751: PetscHashJKKey key;
1752: PetscHashJKIter missing, iter;
1754: key.j = findices[i];
1755: if (key.j >= 0) {
1756: /* Get indices for coarse elements */
1757: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1758: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1759: for (c = 0; c < numCIndices; ++c) {
1760: key.k = cindices[c];
1761: if (key.k < 0) continue;
1762: PetscHashJKPut(ht, key, &missing, &iter);
1763: if (missing) {
1764: PetscHashJKSet(ht, iter, 1);
1765: if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1766: else ++onz[key.j-rStart];
1767: }
1768: }
1769: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1770: }
1771: }
1772: }
1773: PetscSFDestroy(&coarseCellSF);
1774: VecDestroy(&pointVec);
1775: }
1776: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1777: }
1778: }
1779: PetscHashJKDestroy(&ht);
1780: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1781: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1782: PetscFree2(dnz,onz);
1783: for (field = 0; field < Nf; ++field) {
1784: PetscObject obj;
1785: PetscClassId id;
1786: PetscDualSpace Q = NULL;
1787: PetscQuadrature f;
1788: const PetscReal *qpoints, *qweights;
1789: PetscInt Nc, qNc, Np, fpdim, i, d;
1791: PetscDSGetDiscretization(prob, field, &obj);
1792: PetscObjectGetClassId(obj, &id);
1793: if (id == PETSCFE_CLASSID) {
1794: PetscFE fe = (PetscFE) obj;
1796: PetscFEGetDualSpace(fe, &Q);
1797: PetscFEGetNumComponents(fe, &Nc);
1798: } else if (id == PETSCFV_CLASSID) {
1799: PetscFV fv = (PetscFV) obj;
1801: PetscFVGetDualSpace(fv, &Q);
1802: Nc = 1;
1803: } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field);
1804: PetscDualSpaceGetDimension(Q, &fpdim);
1805: /* For each fine grid cell */
1806: for (cell = cStart; cell < cEnd; ++cell) {
1807: PetscInt *findices, *cindices;
1808: PetscInt numFIndices, numCIndices;
1810: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1811: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1812: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1813: for (i = 0; i < fpdim; ++i) {
1814: Vec pointVec;
1815: PetscScalar *pV;
1816: PetscSF coarseCellSF = NULL;
1817: const PetscSFNode *coarseCells;
1818: PetscInt numCoarseCells, cpdim, q, c, j;
1820: /* Get points from the dual basis functional quadrature */
1821: PetscDualSpaceGetFunctional(Q, i, &f);
1822: PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
1823: if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, Nc);
1824: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1825: VecSetBlockSize(pointVec, dim);
1826: VecGetArray(pointVec, &pV);
1827: for (q = 0; q < Np; ++q) {
1828: /* Transform point to real space */
1829: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1830: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1831: }
1832: VecRestoreArray(pointVec, &pV);
1833: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1834: /* OPT: Read this out from preallocation information */
1835: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1836: /* Update preallocation info */
1837: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1838: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1839: VecGetArray(pointVec, &pV);
1840: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1841: PetscReal pVReal[3];
1843: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1844: /* Transform points from real space to coarse reference space */
1845: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
1846: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
1847: CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x);
1849: if (id == PETSCFE_CLASSID) {
1850: PetscFE fe = (PetscFE) obj;
1851: PetscReal *B;
1853: /* Evaluate coarse basis on contained point */
1854: PetscFEGetDimension(fe, &cpdim);
1855: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1856: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
1857: /* Get elemMat entries by multiplying by weight */
1858: for (j = 0; j < cpdim; ++j) {
1859: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
1860: }
1861: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1862: } else {
1863: cpdim = 1;
1864: for (j = 0; j < cpdim; ++j) {
1865: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
1866: }
1867: }
1868: /* Update interpolator */
1869: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
1870: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
1871: MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
1872: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1873: }
1874: VecRestoreArray(pointVec, &pV);
1875: PetscSFDestroy(&coarseCellSF);
1876: VecDestroy(&pointVec);
1877: }
1878: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1879: }
1880: }
1881: PetscFree3(v0,J,invJ);
1882: PetscFree3(v0c,Jc,invJc);
1883: PetscFree(elemMat);
1884: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1885: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1886: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1887: return(0);
1888: }
1890: /*@
1891: DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM.
1893: Input Parameters:
1894: + dmf - The fine mesh
1895: . dmc - The coarse mesh
1896: - user - The user context
1898: Output Parameter:
1899: . mass - The mass matrix
1901: Level: developer
1903: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1904: @*/
1905: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
1906: {
1907: DM_Plex *mesh = (DM_Plex *) dmf->data;
1908: const char *name = "Mass Matrix";
1909: PetscDS prob;
1910: PetscSection fsection, csection, globalFSection, globalCSection;
1911: PetscHashJK ht;
1912: PetscLayout rLayout;
1913: PetscInt *dnz, *onz;
1914: PetscInt locRows, rStart, rEnd;
1915: PetscReal *x, *v0, *J, *invJ, detJ;
1916: PetscReal *v0c, *Jc, *invJc, detJc;
1917: PetscScalar *elemMat;
1918: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
1922: DMGetCoordinateDim(dmc, &dim);
1923: DMGetDS(dmc, &prob);
1924: PetscDSGetRefCoordArrays(prob, &x, NULL);
1925: PetscDSGetNumFields(prob, &Nf);
1926: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1927: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1928: DMGetDefaultSection(dmf, &fsection);
1929: DMGetDefaultGlobalSection(dmf, &globalFSection);
1930: DMGetDefaultSection(dmc, &csection);
1931: DMGetDefaultGlobalSection(dmc, &globalCSection);
1932: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1933: PetscDSGetTotalDimension(prob, &totDim);
1934: PetscMalloc1(totDim, &elemMat);
1936: MatGetLocalSize(mass, &locRows, NULL);
1937: PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
1938: PetscLayoutSetLocalSize(rLayout, locRows);
1939: PetscLayoutSetBlockSize(rLayout, 1);
1940: PetscLayoutSetUp(rLayout);
1941: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1942: PetscLayoutDestroy(&rLayout);
1943: PetscCalloc2(locRows,&dnz,locRows,&onz);
1944: PetscHashJKCreate(&ht);
1945: for (field = 0; field < Nf; ++field) {
1946: PetscObject obj;
1947: PetscClassId id;
1948: PetscQuadrature quad;
1949: const PetscReal *qpoints;
1950: PetscInt Nq, Nc, i, d;
1952: PetscDSGetDiscretization(prob, field, &obj);
1953: PetscObjectGetClassId(obj, &id);
1954: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
1955: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
1956: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
1957: /* For each fine grid cell */
1958: for (cell = cStart; cell < cEnd; ++cell) {
1959: Vec pointVec;
1960: PetscScalar *pV;
1961: PetscSF coarseCellSF = NULL;
1962: const PetscSFNode *coarseCells;
1963: PetscInt numCoarseCells, q, c;
1964: PetscInt *findices, *cindices;
1965: PetscInt numFIndices, numCIndices;
1967: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1968: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1969: /* Get points from the quadrature */
1970: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
1971: VecSetBlockSize(pointVec, dim);
1972: VecGetArray(pointVec, &pV);
1973: for (q = 0; q < Nq; ++q) {
1974: /* Transform point to real space */
1975: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1976: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1977: }
1978: VecRestoreArray(pointVec, &pV);
1979: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1980: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1981: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1982: /* Update preallocation info */
1983: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1984: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1985: {
1986: PetscHashJKKey key;
1987: PetscHashJKIter missing, iter;
1989: for (i = 0; i < numFIndices; ++i) {
1990: key.j = findices[i];
1991: if (key.j >= 0) {
1992: /* Get indices for coarse elements */
1993: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1994: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1995: for (c = 0; c < numCIndices; ++c) {
1996: key.k = cindices[c];
1997: if (key.k < 0) continue;
1998: PetscHashJKPut(ht, key, &missing, &iter);
1999: if (missing) {
2000: PetscHashJKSet(ht, iter, 1);
2001: if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
2002: else ++onz[key.j-rStart];
2003: }
2004: }
2005: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2006: }
2007: }
2008: }
2009: }
2010: PetscSFDestroy(&coarseCellSF);
2011: VecDestroy(&pointVec);
2012: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2013: }
2014: }
2015: PetscHashJKDestroy(&ht);
2016: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2017: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2018: PetscFree2(dnz,onz);
2019: for (field = 0; field < Nf; ++field) {
2020: PetscObject obj;
2021: PetscClassId id;
2022: PetscQuadrature quad;
2023: PetscReal *Bfine;
2024: const PetscReal *qpoints, *qweights;
2025: PetscInt Nq, Nc, i, d;
2027: PetscDSGetDiscretization(prob, field, &obj);
2028: PetscObjectGetClassId(obj, &id);
2029: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);}
2030: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2031: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
2032: /* For each fine grid cell */
2033: for (cell = cStart; cell < cEnd; ++cell) {
2034: Vec pointVec;
2035: PetscScalar *pV;
2036: PetscSF coarseCellSF = NULL;
2037: const PetscSFNode *coarseCells;
2038: PetscInt numCoarseCells, cpdim, q, c, j;
2039: PetscInt *findices, *cindices;
2040: PetscInt numFIndices, numCIndices;
2042: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2043: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2044: /* Get points from the quadrature */
2045: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2046: VecSetBlockSize(pointVec, dim);
2047: VecGetArray(pointVec, &pV);
2048: for (q = 0; q < Nq; ++q) {
2049: /* Transform point to real space */
2050: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
2051: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2052: }
2053: VecRestoreArray(pointVec, &pV);
2054: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2055: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2056: /* Update matrix */
2057: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2058: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2059: VecGetArray(pointVec, &pV);
2060: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2061: PetscReal pVReal[3];
2063: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2064: /* Transform points from real space to coarse reference space */
2065: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2066: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2067: CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x);
2069: if (id == PETSCFE_CLASSID) {
2070: PetscFE fe = (PetscFE) obj;
2071: PetscReal *B;
2073: /* Evaluate coarse basis on contained point */
2074: PetscFEGetDimension(fe, &cpdim);
2075: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2076: /* Get elemMat entries by multiplying by weight */
2077: for (i = 0; i < numFIndices; ++i) {
2078: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2079: for (j = 0; j < cpdim; ++j) {
2080: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2081: }
2082: /* Update interpolator */
2083: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2084: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2085: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2086: }
2087: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2088: } else {
2089: cpdim = 1;
2090: for (i = 0; i < numFIndices; ++i) {
2091: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2092: for (j = 0; j < cpdim; ++j) {
2093: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2094: }
2095: /* Update interpolator */
2096: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2097: PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);
2098: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2099: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2100: }
2101: }
2102: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2103: }
2104: VecRestoreArray(pointVec, &pV);
2105: PetscSFDestroy(&coarseCellSF);
2106: VecDestroy(&pointVec);
2107: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2108: }
2109: }
2110: PetscFree3(v0,J,invJ);
2111: PetscFree3(v0c,Jc,invJc);
2112: PetscFree(elemMat);
2113: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2114: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2115: return(0);
2116: }
2118: /*@
2119: DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns
2121: Input Parameters:
2122: + dmc - The coarse mesh
2123: - dmf - The fine mesh
2124: - user - The user context
2126: Output Parameter:
2127: . sc - The mapping
2129: Level: developer
2131: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2132: @*/
2133: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
2134: {
2135: PetscDS prob;
2136: PetscFE *feRef;
2137: PetscFV *fvRef;
2138: Vec fv, cv;
2139: IS fis, cis;
2140: PetscSection fsection, fglobalSection, csection, cglobalSection;
2141: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
2142: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
2146: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2147: DMGetDimension(dmf, &dim);
2148: DMGetDefaultSection(dmf, &fsection);
2149: DMGetDefaultGlobalSection(dmf, &fglobalSection);
2150: DMGetDefaultSection(dmc, &csection);
2151: DMGetDefaultGlobalSection(dmc, &cglobalSection);
2152: PetscSectionGetNumFields(fsection, &Nf);
2153: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2154: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2155: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2156: DMGetDS(dmc, &prob);
2157: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
2158: for (f = 0; f < Nf; ++f) {
2159: PetscObject obj;
2160: PetscClassId id;
2161: PetscInt fNb = 0, Nc = 0;
2163: PetscDSGetDiscretization(prob, f, &obj);
2164: PetscObjectGetClassId(obj, &id);
2165: if (id == PETSCFE_CLASSID) {
2166: PetscFE fe = (PetscFE) obj;
2168: PetscFERefine(fe, &feRef[f]);
2169: PetscFEGetDimension(feRef[f], &fNb);
2170: PetscFEGetNumComponents(fe, &Nc);
2171: } else if (id == PETSCFV_CLASSID) {
2172: PetscFV fv = (PetscFV) obj;
2173: PetscDualSpace Q;
2175: PetscFVRefine(fv, &fvRef[f]);
2176: PetscFVGetDualSpace(fvRef[f], &Q);
2177: PetscDualSpaceGetDimension(Q, &fNb);
2178: PetscFVGetNumComponents(fv, &Nc);
2179: }
2180: fTotDim += fNb*Nc;
2181: }
2182: PetscDSGetTotalDimension(prob, &cTotDim);
2183: PetscMalloc1(cTotDim,&cmap);
2184: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
2185: PetscFE feC;
2186: PetscFV fvC;
2187: PetscDualSpace QF, QC;
2188: PetscInt NcF, NcC, fpdim, cpdim;
2190: if (feRef[field]) {
2191: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
2192: PetscFEGetNumComponents(feC, &NcC);
2193: PetscFEGetNumComponents(feRef[field], &NcF);
2194: PetscFEGetDualSpace(feRef[field], &QF);
2195: PetscDualSpaceGetDimension(QF, &fpdim);
2196: PetscFEGetDualSpace(feC, &QC);
2197: PetscDualSpaceGetDimension(QC, &cpdim);
2198: } else {
2199: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
2200: PetscFVGetNumComponents(fvC, &NcC);
2201: PetscFVGetNumComponents(fvRef[field], &NcF);
2202: PetscFVGetDualSpace(fvRef[field], &QF);
2203: PetscDualSpaceGetDimension(QF, &fpdim);
2204: PetscFVGetDualSpace(fvC, &QC);
2205: PetscDualSpaceGetDimension(QC, &cpdim);
2206: }
2207: 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);
2208: for (c = 0; c < cpdim; ++c) {
2209: PetscQuadrature cfunc;
2210: const PetscReal *cqpoints;
2211: PetscInt NpC;
2212: PetscBool found = PETSC_FALSE;
2214: PetscDualSpaceGetFunctional(QC, c, &cfunc);
2215: PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);
2216: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
2217: for (f = 0; f < fpdim; ++f) {
2218: PetscQuadrature ffunc;
2219: const PetscReal *fqpoints;
2220: PetscReal sum = 0.0;
2221: PetscInt NpF, comp;
2223: PetscDualSpaceGetFunctional(QF, f, &ffunc);
2224: PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);
2225: if (NpC != NpF) continue;
2226: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
2227: if (sum > 1.0e-9) continue;
2228: for (comp = 0; comp < NcC; ++comp) {
2229: cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
2230: }
2231: found = PETSC_TRUE;
2232: break;
2233: }
2234: if (!found) {
2235: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
2236: if (fvRef[field]) {
2237: PetscInt comp;
2238: for (comp = 0; comp < NcC; ++comp) {
2239: cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
2240: }
2241: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
2242: }
2243: }
2244: offsetC += cpdim*NcC;
2245: offsetF += fpdim*NcF;
2246: }
2247: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
2248: PetscFree2(feRef,fvRef);
2250: DMGetGlobalVector(dmf, &fv);
2251: DMGetGlobalVector(dmc, &cv);
2252: VecGetOwnershipRange(cv, &startC, &endC);
2253: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
2254: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
2255: PetscMalloc1(m,&cindices);
2256: PetscMalloc1(m,&findices);
2257: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
2258: for (c = cStart; c < cEnd; ++c) {
2259: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
2260: for (d = 0; d < cTotDim; ++d) {
2261: if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
2262: 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]]);
2263: cindices[cellCIndices[d]-startC] = cellCIndices[d];
2264: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
2265: }
2266: }
2267: PetscFree(cmap);
2268: PetscFree2(cellCIndices,cellFIndices);
2270: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
2271: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
2272: VecScatterCreate(cv, cis, fv, fis, sc);
2273: ISDestroy(&cis);
2274: ISDestroy(&fis);
2275: DMRestoreGlobalVector(dmf, &fv);
2276: DMRestoreGlobalVector(dmc, &cv);
2277: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2278: return(0);
2279: }