Actual source code: plexfem.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
4: #include <petsc/private/hashsetij.h>
5: #include <petsc/private/petscfeimpl.h>
6: #include <petsc/private/petscfvimpl.h>
8: static PetscErrorCode DMPlexConvertPlex(DM dm, DM *plex, PetscBool copy)
9: {
10: PetscBool isPlex;
14: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
15: if (isPlex) {
16: *plex = dm;
17: PetscObjectReference((PetscObject) dm);
18: } else {
19: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
20: if (!*plex) {
21: DMConvert(dm,DMPLEX,plex);
22: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
23: if (copy) {
24: const char *comps[3] = {"A", "dmAux"};
25: PetscObject obj;
26: PetscInt i;
28: {
29: /* Run the subdomain hook (this will copy the DMSNES/DMTS) */
30: DMSubDomainHookLink link;
31: for (link = dm->subdomainhook; link; link = link->next) {
32: if (link->ddhook) {(*link->ddhook)(dm, *plex, link->ctx);}
33: }
34: }
35: for (i = 0; i < 3; i++) {
36: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
37: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
38: }
39: }
40: } else {
41: PetscObjectReference((PetscObject) *plex);
42: }
43: }
44: return(0);
45: }
47: static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom (void *ctx)
48: {
49: PetscFEGeom *geom = (PetscFEGeom *) ctx;
53: PetscFEGeomDestroy(&geom);
54: return(0);
55: }
57: static PetscErrorCode DMPlexGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
58: {
59: char composeStr[33] = {0};
60: PetscObjectId id;
61: PetscContainer container;
62: PetscErrorCode ierr;
65: PetscObjectGetId((PetscObject)quad,&id);
66: PetscSNPrintf(composeStr, 32, "DMPlexGetFEGeom_%x\n", id);
67: PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
68: if (container) {
69: PetscContainerGetPointer(container, (void **) geom);
70: } else {
71: DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
72: PetscContainerCreate(PETSC_COMM_SELF,&container);
73: PetscContainerSetPointer(container, (void *) *geom);
74: PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
75: PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
76: PetscContainerDestroy(&container);
77: }
78: return(0);
79: }
81: static PetscErrorCode DMPlexRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
82: {
84: *geom = NULL;
85: return(0);
86: }
88: /*@
89: DMPlexGetScale - Get the scale for the specified fundamental unit
91: Not collective
93: Input Arguments:
94: + dm - the DM
95: - unit - The SI unit
97: Output Argument:
98: . scale - The value used to scale all quantities with this unit
100: Level: advanced
102: .seealso: DMPlexSetScale(), PetscUnit
103: @*/
104: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
105: {
106: DM_Plex *mesh = (DM_Plex*) dm->data;
111: *scale = mesh->scale[unit];
112: return(0);
113: }
115: /*@
116: DMPlexSetScale - Set the scale for the specified fundamental unit
118: Not collective
120: Input Arguments:
121: + dm - the DM
122: . unit - The SI unit
123: - scale - The value used to scale all quantities with this unit
125: Level: advanced
127: .seealso: DMPlexGetScale(), PetscUnit
128: @*/
129: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
130: {
131: DM_Plex *mesh = (DM_Plex*) dm->data;
135: mesh->scale[unit] = scale;
136: return(0);
137: }
139: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
140: {
141: 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}}};
142: PetscInt *ctxInt = (PetscInt *) ctx;
143: PetscInt dim2 = ctxInt[0];
144: PetscInt d = ctxInt[1];
145: PetscInt i, j, k = dim > 2 ? d - dim : d;
148: if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
149: for (i = 0; i < dim; i++) mode[i] = 0.;
150: if (d < dim) {
151: mode[d] = 1.; /* Translation along axis d */
152: } else {
153: for (i = 0; i < dim; i++) {
154: for (j = 0; j < dim; j++) {
155: mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */
156: }
157: }
158: }
159: return(0);
160: }
162: /*@
163: DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation
165: Collective on DM
167: Input Arguments:
168: . dm - the DM
170: Output Argument:
171: . sp - the null space
173: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
175: Level: advanced
177: .seealso: MatNullSpaceCreate(), PCGAMG
178: @*/
179: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
180: {
181: MPI_Comm comm;
182: Vec mode[6];
183: PetscSection section, globalSection;
184: PetscInt dim, dimEmbed, n, m, mmin, d, i, j;
188: PetscObjectGetComm((PetscObject)dm,&comm);
189: DMGetDimension(dm, &dim);
190: DMGetCoordinateDim(dm, &dimEmbed);
191: if (dim == 1) {
192: MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
193: return(0);
194: }
195: DMGetSection(dm, §ion);
196: DMGetGlobalSection(dm, &globalSection);
197: PetscSectionGetConstrainedStorageSize(globalSection, &n);
198: m = (dim*(dim+1))/2;
199: VecCreate(comm, &mode[0]);
200: VecSetSizes(mode[0], n, PETSC_DETERMINE);
201: VecSetUp(mode[0]);
202: VecGetSize(mode[0], &n);
203: mmin = PetscMin(m, n);
204: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
205: for (d = 0; d < m; d++) {
206: PetscInt ctx[2];
207: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
208: void *voidctx = (void *) (&ctx[0]);
210: ctx[0] = dimEmbed;
211: ctx[1] = d;
212: DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
213: }
214: for (i = 0; i < PetscMin(dim, mmin); ++i) {VecNormalize(mode[i], NULL);}
215: /* Orthonormalize system */
216: for (i = dim; i < mmin; ++i) {
217: PetscScalar dots[6];
219: VecMDot(mode[i], i, mode, dots);
220: for (j = 0; j < i; ++j) dots[j] *= -1.0;
221: VecMAXPY(mode[i], i, dots, mode);
222: VecNormalize(mode[i], NULL);
223: }
224: MatNullSpaceCreate(comm, PETSC_FALSE, mmin, mode, sp);
225: for (i = 0; i < m; ++i) {VecDestroy(&mode[i]);}
226: return(0);
227: }
229: /*@
230: DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation
232: Collective on DM
234: Input Arguments:
235: + dm - the DM
236: . nb - The number of bodies
237: . label - The DMLabel marking each domain
238: . nids - The number of ids per body
239: - ids - An array of the label ids in sequence for each domain
241: Output Argument:
242: . sp - the null space
244: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
246: Level: advanced
248: .seealso: MatNullSpaceCreate()
249: @*/
250: PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp)
251: {
252: MPI_Comm comm;
253: PetscSection section, globalSection;
254: Vec *mode;
255: PetscScalar *dots;
256: PetscInt dim, dimEmbed, n, m, b, d, i, j, off;
260: PetscObjectGetComm((PetscObject)dm,&comm);
261: DMGetDimension(dm, &dim);
262: DMGetCoordinateDim(dm, &dimEmbed);
263: DMGetSection(dm, §ion);
264: DMGetGlobalSection(dm, &globalSection);
265: PetscSectionGetConstrainedStorageSize(globalSection, &n);
266: m = nb * (dim*(dim+1))/2;
267: PetscMalloc2(m, &mode, m, &dots);
268: VecCreate(comm, &mode[0]);
269: VecSetSizes(mode[0], n, PETSC_DETERMINE);
270: VecSetUp(mode[0]);
271: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
272: for (b = 0, off = 0; b < nb; ++b) {
273: for (d = 0; d < m/nb; ++d) {
274: PetscInt ctx[2];
275: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
276: void *voidctx = (void *) (&ctx[0]);
278: ctx[0] = dimEmbed;
279: ctx[1] = d;
280: DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
281: off += nids[b];
282: }
283: }
284: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
285: /* Orthonormalize system */
286: for (i = 0; i < m; ++i) {
287: VecMDot(mode[i], i, mode, dots);
288: for (j = 0; j < i; ++j) dots[j] *= -1.0;
289: VecMAXPY(mode[i], i, dots, mode);
290: VecNormalize(mode[i], NULL);
291: }
292: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
293: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
294: PetscFree2(mode, dots);
295: return(0);
296: }
298: /*@
299: DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
300: are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
301: evaluating the dual space basis of that point. A basis function is associated with the point in its
302: transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
303: projection height, which is set with this function. By default, the maximum projection height is zero, which means
304: that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior
305: basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
307: Input Parameters:
308: + dm - the DMPlex object
309: - height - the maximum projection height >= 0
311: Level: advanced
313: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
314: @*/
315: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
316: {
317: DM_Plex *plex = (DM_Plex *) dm->data;
321: plex->maxProjectionHeight = height;
322: return(0);
323: }
325: /*@
326: DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
327: DMPlexProjectXXXLocal() functions.
329: Input Parameters:
330: . dm - the DMPlex object
332: Output Parameters:
333: . height - the maximum projection height
335: Level: intermediate
337: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
338: @*/
339: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
340: {
341: DM_Plex *plex = (DM_Plex *) dm->data;
345: *height = plex->maxProjectionHeight;
346: return(0);
347: }
349: /*@C
350: DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector
352: Input Parameters:
353: + dm - The DM, with a PetscDS that matches the problem being constrained
354: . time - The time
355: . field - The field to constrain
356: . Nc - The number of constrained field components, or 0 for all components
357: . comps - An array of constrained component numbers, or NULL for all components
358: . label - The DMLabel defining constrained points
359: . numids - The number of DMLabel ids for constrained points
360: . ids - An array of ids for constrained points
361: . func - A pointwise function giving boundary values
362: - ctx - An optional user context for bcFunc
364: Output Parameter:
365: . locX - A local vector to receives the boundary values
367: Level: developer
369: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
370: @*/
371: 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)
372: {
373: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
374: void **ctxs;
375: PetscInt numFields;
376: PetscErrorCode ierr;
379: DMGetNumFields(dm, &numFields);
380: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
381: funcs[field] = func;
382: ctxs[field] = ctx;
383: DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
384: PetscFree2(funcs,ctxs);
385: return(0);
386: }
388: /*@C
389: DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector
391: Input Parameters:
392: + dm - The DM, with a PetscDS that matches the problem being constrained
393: . time - The time
394: . locU - A local vector with the input solution values
395: . field - The field to constrain
396: . Nc - The number of constrained field components, or 0 for all components
397: . comps - An array of constrained component numbers, or NULL for all components
398: . label - The DMLabel defining constrained points
399: . numids - The number of DMLabel ids for constrained points
400: . ids - An array of ids for constrained points
401: . func - A pointwise function giving boundary values
402: - ctx - An optional user context for bcFunc
404: Output Parameter:
405: . locX - A local vector to receives the boundary values
407: Level: developer
409: .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary()
410: @*/
411: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
412: void (*func)(PetscInt, PetscInt, PetscInt,
413: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
414: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
415: PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
416: PetscScalar[]),
417: void *ctx, Vec locX)
418: {
419: void (**funcs)(PetscInt, PetscInt, PetscInt,
420: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
421: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
422: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
423: void **ctxs;
424: PetscInt numFields;
425: PetscErrorCode ierr;
428: DMGetNumFields(dm, &numFields);
429: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
430: funcs[field] = func;
431: ctxs[field] = ctx;
432: DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
433: PetscFree2(funcs,ctxs);
434: return(0);
435: }
437: /*@C
438: DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector
440: Input Parameters:
441: + dm - The DM, with a PetscDS that matches the problem being constrained
442: . time - The time
443: . faceGeometry - A vector with the FVM face geometry information
444: . cellGeometry - A vector with the FVM cell geometry information
445: . Grad - A vector with the FVM cell gradient information
446: . field - The field to constrain
447: . Nc - The number of constrained field components, or 0 for all components
448: . comps - An array of constrained component numbers, or NULL for all components
449: . label - The DMLabel defining constrained points
450: . numids - The number of DMLabel ids for constrained points
451: . ids - An array of ids for constrained points
452: . func - A pointwise function giving boundary values
453: - ctx - An optional user context for bcFunc
455: Output Parameter:
456: . locX - A local vector to receives the boundary values
458: Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary()
460: Level: developer
462: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
463: @*/
464: 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[],
465: PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
466: {
467: PetscDS prob;
468: PetscSF sf;
469: DM dmFace, dmCell, dmGrad;
470: const PetscScalar *facegeom, *cellgeom = NULL, *grad;
471: const PetscInt *leaves;
472: PetscScalar *x, *fx;
473: PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i;
474: PetscErrorCode ierr, ierru = 0;
477: DMGetPointSF(dm, &sf);
478: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
479: nleaves = PetscMax(0, nleaves);
480: DMGetDimension(dm, &dim);
481: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
482: DMGetDS(dm, &prob);
483: VecGetDM(faceGeometry, &dmFace);
484: VecGetArrayRead(faceGeometry, &facegeom);
485: if (cellGeometry) {
486: VecGetDM(cellGeometry, &dmCell);
487: VecGetArrayRead(cellGeometry, &cellgeom);
488: }
489: if (Grad) {
490: PetscFV fv;
492: PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
493: VecGetDM(Grad, &dmGrad);
494: VecGetArrayRead(Grad, &grad);
495: PetscFVGetNumComponents(fv, &pdim);
496: DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
497: }
498: VecGetArray(locX, &x);
499: for (i = 0; i < numids; ++i) {
500: IS faceIS;
501: const PetscInt *faces;
502: PetscInt numFaces, f;
504: DMLabelGetStratumIS(label, ids[i], &faceIS);
505: if (!faceIS) continue; /* No points with that id on this process */
506: ISGetLocalSize(faceIS, &numFaces);
507: ISGetIndices(faceIS, &faces);
508: for (f = 0; f < numFaces; ++f) {
509: const PetscInt face = faces[f], *cells;
510: PetscFVFaceGeom *fg;
512: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
513: PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
514: if (loc >= 0) continue;
515: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
516: DMPlexGetSupport(dm, face, &cells);
517: if (Grad) {
518: PetscFVCellGeom *cg;
519: PetscScalar *cx, *cgrad;
520: PetscScalar *xG;
521: PetscReal dx[3];
522: PetscInt d;
524: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
525: DMPlexPointLocalRead(dm, cells[0], x, &cx);
526: DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
527: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
528: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
529: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
530: ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
531: if (ierru) {
532: ISRestoreIndices(faceIS, &faces);
533: ISDestroy(&faceIS);
534: goto cleanup;
535: }
536: } else {
537: PetscScalar *xI;
538: PetscScalar *xG;
540: DMPlexPointLocalRead(dm, cells[0], x, &xI);
541: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
542: ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
543: if (ierru) {
544: ISRestoreIndices(faceIS, &faces);
545: ISDestroy(&faceIS);
546: goto cleanup;
547: }
548: }
549: }
550: ISRestoreIndices(faceIS, &faces);
551: ISDestroy(&faceIS);
552: }
553: cleanup:
554: VecRestoreArray(locX, &x);
555: if (Grad) {
556: DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
557: VecRestoreArrayRead(Grad, &grad);
558: }
559: if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
560: VecRestoreArrayRead(faceGeometry, &facegeom);
561: CHKERRQ(ierru);
562: return(0);
563: }
565: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
566: {
567: PetscDS prob;
568: PetscInt numBd, b;
572: DMGetDS(dm, &prob);
573: PetscDSGetNumBoundary(prob, &numBd);
574: for (b = 0; b < numBd; ++b) {
575: DMBoundaryConditionType type;
576: const char *labelname;
577: DMLabel label;
578: PetscInt field, Nc;
579: const PetscInt *comps;
580: PetscObject obj;
581: PetscClassId id;
582: void (*func)(void);
583: PetscInt numids;
584: const PetscInt *ids;
585: void *ctx;
587: DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
588: if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
589: DMGetLabel(dm, labelname, &label);
590: DMGetField(dm, field, NULL, &obj);
591: PetscObjectGetClassId(obj, &id);
592: if (id == PETSCFE_CLASSID) {
593: switch (type) {
594: /* for FEM, there is no insertion to be done for non-essential boundary conditions */
595: case DM_BC_ESSENTIAL:
596: DMPlexLabelAddCells(dm,label);
597: DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
598: DMPlexLabelClearCells(dm,label);
599: break;
600: case DM_BC_ESSENTIAL_FIELD:
601: DMPlexLabelAddCells(dm,label);
602: DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
603: (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
604: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
605: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
606: DMPlexLabelClearCells(dm,label);
607: break;
608: default: break;
609: }
610: } else if (id == PETSCFV_CLASSID) {
611: if (!faceGeomFVM) continue;
612: DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
613: (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
614: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
615: }
616: return(0);
617: }
619: /*@
620: DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector
622: Input Parameters:
623: + dm - The DM
624: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
625: . time - The time
626: . faceGeomFVM - Face geometry data for FV discretizations
627: . cellGeomFVM - Cell geometry data for FV discretizations
628: - gradFVM - Gradient reconstruction data for FV discretizations
630: Output Parameters:
631: . locX - Solution updated with boundary values
633: Level: developer
635: .seealso: DMProjectFunctionLabelLocal()
636: @*/
637: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
638: {
647: PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
648: return(0);
649: }
651: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
652: {
653: Vec localX;
654: PetscErrorCode ierr;
657: DMGetLocalVector(dm, &localX);
658: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
659: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
660: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
661: DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
662: DMRestoreLocalVector(dm, &localX);
663: return(0);
664: }
666: /*@C
667: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
669: Input Parameters:
670: + dm - The DM
671: . time - The time
672: . funcs - The functions to evaluate for each field component
673: . ctxs - Optional array of contexts to pass to each function, or NULL.
674: - localX - The coefficient vector u_h, a local vector
676: Output Parameter:
677: . diff - The diff ||u - u_h||_2
679: Level: developer
681: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
682: @*/
683: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
684: {
685: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
686: PetscSection section;
687: PetscQuadrature quad;
688: PetscScalar *funcVal, *interpolant;
689: PetscReal *coords, *detJ, *J;
690: PetscReal localDiff = 0.0;
691: const PetscReal *quadWeights;
692: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset;
693: PetscErrorCode ierr;
696: DMGetDimension(dm, &dim);
697: DMGetCoordinateDim(dm, &coordDim);
698: DMGetSection(dm, §ion);
699: PetscSectionGetNumFields(section, &numFields);
700: for (field = 0; field < numFields; ++field) {
701: PetscObject obj;
702: PetscClassId id;
703: PetscInt Nc;
705: DMGetField(dm, field, NULL, &obj);
706: PetscObjectGetClassId(obj, &id);
707: if (id == PETSCFE_CLASSID) {
708: PetscFE fe = (PetscFE) obj;
710: PetscFEGetQuadrature(fe, &quad);
711: PetscFEGetNumComponents(fe, &Nc);
712: } else if (id == PETSCFV_CLASSID) {
713: PetscFV fv = (PetscFV) obj;
715: PetscFVGetQuadrature(fv, &quad);
716: PetscFVGetNumComponents(fv, &Nc);
717: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
718: numComponents += Nc;
719: }
720: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
721: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
722: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
723: DMPlexGetVTKCellHeight(dm, &cellHeight);
724: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
725: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
726: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
727: for (c = cStart; c < cEnd; ++c) {
728: PetscScalar *x = NULL;
729: PetscReal elemDiff = 0.0;
730: PetscInt qc = 0;
732: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
733: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
735: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
736: PetscObject obj;
737: PetscClassId id;
738: void * const ctx = ctxs ? ctxs[field] : NULL;
739: PetscInt Nb, Nc, q, fc;
741: DMGetField(dm, field, NULL, &obj);
742: PetscObjectGetClassId(obj, &id);
743: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
744: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
745: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
746: if (debug) {
747: char title[1024];
748: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
749: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
750: }
751: for (q = 0; q < Nq; ++q) {
752: 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);
753: (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
754: if (ierr) {
755: PetscErrorCode ierr2;
756: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
757: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
758: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
759:
760: }
761: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
762: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
763: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
764: for (fc = 0; fc < Nc; ++fc) {
765: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
766: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
767: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
768: }
769: }
770: fieldOffset += Nb;
771: qc += Nc;
772: }
773: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
774: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
775: localDiff += elemDiff;
776: }
777: PetscFree5(funcVal,interpolant,coords,detJ,J);
778: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
779: *diff = PetscSqrtReal(*diff);
780: return(0);
781: }
783: 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)
784: {
785: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
786: PetscSection section;
787: PetscQuadrature quad;
788: Vec localX;
789: PetscScalar *funcVal, *interpolant;
790: const PetscReal *quadPoints, *quadWeights;
791: PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ;
792: PetscReal localDiff = 0.0;
793: PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
794: PetscErrorCode ierr;
797: DMGetDimension(dm, &dim);
798: DMGetCoordinateDim(dm, &coordDim);
799: DMGetSection(dm, §ion);
800: PetscSectionGetNumFields(section, &numFields);
801: DMGetLocalVector(dm, &localX);
802: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
803: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
804: for (field = 0; field < numFields; ++field) {
805: PetscFE fe;
806: PetscInt Nc;
808: DMGetField(dm, field, NULL, (PetscObject *) &fe);
809: PetscFEGetQuadrature(fe, &quad);
810: PetscFEGetNumComponents(fe, &Nc);
811: numComponents += Nc;
812: }
813: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
814: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
815: /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
816: PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
817: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
818: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
819: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
820: for (c = cStart; c < cEnd; ++c) {
821: PetscScalar *x = NULL;
822: PetscReal elemDiff = 0.0;
823: PetscInt qc = 0;
825: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);
826: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
828: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
829: PetscFE fe;
830: void * const ctx = ctxs ? ctxs[field] : NULL;
831: PetscReal *basisDer;
832: PetscInt Nb, Nc, q, fc;
834: DMGetField(dm, field, NULL, (PetscObject *) &fe);
835: PetscFEGetDimension(fe, &Nb);
836: PetscFEGetNumComponents(fe, &Nc);
837: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
838: if (debug) {
839: char title[1024];
840: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
841: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
842: }
843: for (q = 0; q < Nq; ++q) {
844: 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);
845: (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
846: if (ierr) {
847: PetscErrorCode ierr2;
848: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
849: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
850: ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
851:
852: }
853: PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
854: for (fc = 0; fc < Nc; ++fc) {
855: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
856: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
857: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
858: }
859: }
860: fieldOffset += Nb;
861: qc += Nc;
862: }
863: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
864: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
865: localDiff += elemDiff;
866: }
867: PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
868: DMRestoreLocalVector(dm, &localX);
869: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
870: *diff = PetscSqrtReal(*diff);
871: return(0);
872: }
874: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
875: {
876: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
877: PetscSection section;
878: PetscQuadrature quad;
879: Vec localX;
880: PetscScalar *funcVal, *interpolant;
881: PetscReal *coords, *detJ, *J;
882: PetscReal *localDiff;
883: const PetscReal *quadPoints, *quadWeights;
884: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
885: PetscErrorCode ierr;
888: DMGetDimension(dm, &dim);
889: DMGetCoordinateDim(dm, &coordDim);
890: DMGetSection(dm, §ion);
891: PetscSectionGetNumFields(section, &numFields);
892: DMGetLocalVector(dm, &localX);
893: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
894: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
895: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
896: for (field = 0; field < numFields; ++field) {
897: PetscObject obj;
898: PetscClassId id;
899: PetscInt Nc;
901: DMGetField(dm, field, NULL, &obj);
902: PetscObjectGetClassId(obj, &id);
903: if (id == PETSCFE_CLASSID) {
904: PetscFE fe = (PetscFE) obj;
906: PetscFEGetQuadrature(fe, &quad);
907: PetscFEGetNumComponents(fe, &Nc);
908: } else if (id == PETSCFV_CLASSID) {
909: PetscFV fv = (PetscFV) obj;
911: PetscFVGetQuadrature(fv, &quad);
912: PetscFVGetNumComponents(fv, &Nc);
913: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
914: numComponents += Nc;
915: }
916: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
917: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
918: PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
919: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
920: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
921: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
922: for (c = cStart; c < cEnd; ++c) {
923: PetscScalar *x = NULL;
924: PetscInt qc = 0;
926: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
927: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
929: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
930: PetscObject obj;
931: PetscClassId id;
932: void * const ctx = ctxs ? ctxs[field] : NULL;
933: PetscInt Nb, Nc, q, fc;
935: PetscReal elemDiff = 0.0;
937: DMGetField(dm, field, NULL, &obj);
938: PetscObjectGetClassId(obj, &id);
939: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
940: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
941: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
942: if (debug) {
943: char title[1024];
944: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
945: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
946: }
947: for (q = 0; q < Nq; ++q) {
948: 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);
949: (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
950: if (ierr) {
951: PetscErrorCode ierr2;
952: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
953: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
954: ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
955:
956: }
957: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
958: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
959: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
960: for (fc = 0; fc < Nc; ++fc) {
961: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
962: 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]);}
963: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
964: }
965: }
966: fieldOffset += Nb;
967: qc += Nc;
968: localDiff[field] += elemDiff;
969: }
970: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
971: }
972: DMRestoreLocalVector(dm, &localX);
973: MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
974: for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
975: PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
976: return(0);
977: }
979: /*@C
980: 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.
982: Input Parameters:
983: + dm - The DM
984: . time - The time
985: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
986: . ctxs - Optional array of contexts to pass to each function, or NULL.
987: - X - The coefficient vector u_h
989: Output Parameter:
990: . D - A Vec which holds the difference ||u - u_h||_2 for each cell
992: Level: developer
994: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
995: @*/
996: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
997: {
998: PetscSection section;
999: PetscQuadrature quad;
1000: Vec localX;
1001: PetscScalar *funcVal, *interpolant;
1002: PetscReal *coords, *detJ, *J;
1003: const PetscReal *quadPoints, *quadWeights;
1004: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1005: PetscErrorCode ierr;
1008: VecSet(D, 0.0);
1009: DMGetDimension(dm, &dim);
1010: DMGetCoordinateDim(dm, &coordDim);
1011: DMGetSection(dm, §ion);
1012: PetscSectionGetNumFields(section, &numFields);
1013: DMGetLocalVector(dm, &localX);
1014: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1015: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1016: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1017: for (field = 0; field < numFields; ++field) {
1018: PetscObject obj;
1019: PetscClassId id;
1020: PetscInt Nc;
1022: DMGetField(dm, field, NULL, &obj);
1023: PetscObjectGetClassId(obj, &id);
1024: if (id == PETSCFE_CLASSID) {
1025: PetscFE fe = (PetscFE) obj;
1027: PetscFEGetQuadrature(fe, &quad);
1028: PetscFEGetNumComponents(fe, &Nc);
1029: } else if (id == PETSCFV_CLASSID) {
1030: PetscFV fv = (PetscFV) obj;
1032: PetscFVGetQuadrature(fv, &quad);
1033: PetscFVGetNumComponents(fv, &Nc);
1034: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1035: numComponents += Nc;
1036: }
1037: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1038: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1039: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
1040: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1041: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1042: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1043: for (c = cStart; c < cEnd; ++c) {
1044: PetscScalar *x = NULL;
1045: PetscScalar elemDiff = 0.0;
1046: PetscInt qc = 0;
1048: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
1049: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1051: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1052: PetscObject obj;
1053: PetscClassId id;
1054: void * const ctx = ctxs ? ctxs[field] : NULL;
1055: PetscInt Nb, Nc, q, fc;
1057: DMGetField(dm, field, NULL, &obj);
1058: PetscObjectGetClassId(obj, &id);
1059: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1060: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1061: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1062: if (funcs[field]) {
1063: for (q = 0; q < Nq; ++q) {
1064: 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);
1065: (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
1066: if (ierr) {
1067: PetscErrorCode ierr2;
1068: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1069: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
1070: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1071:
1072: }
1073: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1074: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1075: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1076: for (fc = 0; fc < Nc; ++fc) {
1077: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1078: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
1079: }
1080: }
1081: }
1082: fieldOffset += Nb;
1083: qc += Nc;
1084: }
1085: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1086: VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1087: }
1088: PetscFree5(funcVal,interpolant,coords,detJ,J);
1089: DMRestoreLocalVector(dm, &localX);
1090: VecSqrtAbs(D);
1091: return(0);
1092: }
1094: /*@C
1095: DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec.
1097: Input Parameters:
1098: + dm - The DM
1099: - LocX - The coefficient vector u_h
1101: Output Parameter:
1102: . locC - A Vec which holds the Clement interpolant of the gradient
1104: Notes:
1105: Add citation to (Clement, 1975) and definition of the interpolant
1106: \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
1108: Level: developer
1110: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1111: @*/
1112: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1113: {
1114: DM_Plex *mesh = (DM_Plex *) dm->data;
1115: PetscInt debug = mesh->printFEM;
1116: DM dmC;
1117: PetscSection section;
1118: PetscQuadrature quad;
1119: PetscScalar *interpolant, *gradsum;
1120: PetscReal *coords, *detJ, *J, *invJ;
1121: const PetscReal *quadPoints, *quadWeights;
1122: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset;
1123: PetscErrorCode ierr;
1126: VecGetDM(locC, &dmC);
1127: VecSet(locC, 0.0);
1128: DMGetDimension(dm, &dim);
1129: DMGetCoordinateDim(dm, &coordDim);
1130: DMGetSection(dm, §ion);
1131: PetscSectionGetNumFields(section, &numFields);
1132: for (field = 0; field < numFields; ++field) {
1133: PetscObject obj;
1134: PetscClassId id;
1135: PetscInt Nc;
1137: DMGetField(dm, field, NULL, &obj);
1138: PetscObjectGetClassId(obj, &id);
1139: if (id == PETSCFE_CLASSID) {
1140: PetscFE fe = (PetscFE) obj;
1142: PetscFEGetQuadrature(fe, &quad);
1143: PetscFEGetNumComponents(fe, &Nc);
1144: } else if (id == PETSCFV_CLASSID) {
1145: PetscFV fv = (PetscFV) obj;
1147: PetscFVGetQuadrature(fv, &quad);
1148: PetscFVGetNumComponents(fv, &Nc);
1149: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1150: numComponents += Nc;
1151: }
1152: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1153: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1154: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);
1155: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1156: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1157: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1158: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1159: for (v = vStart; v < vEnd; ++v) {
1160: PetscScalar volsum = 0.0;
1161: PetscInt *star = NULL;
1162: PetscInt starSize, st, d, fc;
1164: PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));
1165: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1166: for (st = 0; st < starSize*2; st += 2) {
1167: const PetscInt cell = star[st];
1168: PetscScalar *grad = &gradsum[coordDim*numComponents];
1169: PetscScalar *x = NULL;
1170: PetscReal vol = 0.0;
1172: if ((cell < cStart) || (cell >= cEnd)) continue;
1173: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);
1174: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1175: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1176: PetscObject obj;
1177: PetscClassId id;
1178: PetscInt Nb, Nc, q, qc = 0;
1180: PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));
1181: DMGetField(dm, field, NULL, &obj);
1182: PetscObjectGetClassId(obj, &id);
1183: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1184: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1185: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1186: for (q = 0; q < Nq; ++q) {
1187: 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);
1188: if (ierr) {
1189: PetscErrorCode ierr2;
1190: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1191: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1192: ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2);
1193:
1194: }
1195: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);}
1196: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1197: for (fc = 0; fc < Nc; ++fc) {
1198: const PetscReal wt = quadWeights[q*qNc+qc+fc];
1200: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q];
1201: }
1202: vol += quadWeights[q*qNc]*detJ[q];
1203: }
1204: fieldOffset += Nb;
1205: qc += Nc;
1206: }
1207: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1208: for (fc = 0; fc < numComponents; ++fc) {
1209: for (d = 0; d < coordDim; ++d) {
1210: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1211: }
1212: }
1213: volsum += vol;
1214: if (debug) {
1215: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1216: for (fc = 0; fc < numComponents; ++fc) {
1217: for (d = 0; d < coordDim; ++d) {
1218: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1219: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1220: }
1221: }
1222: PetscPrintf(PETSC_COMM_SELF, "]\n");
1223: }
1224: }
1225: for (fc = 0; fc < numComponents; ++fc) {
1226: for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1227: }
1228: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1229: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1230: }
1231: PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);
1232: return(0);
1233: }
1235: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1236: {
1237: DM dmAux = NULL;
1238: PetscDS prob, probAux = NULL;
1239: PetscSection section, sectionAux;
1240: Vec locX, locA;
1241: PetscInt dim, numCells = cEnd - cStart, c, f;
1242: PetscBool useFVM = PETSC_FALSE;
1243: /* DS */
1244: PetscInt Nf, totDim, *uOff, *uOff_x, numConstants;
1245: PetscInt NfAux, totDimAux, *aOff;
1246: PetscScalar *u, *a;
1247: const PetscScalar *constants;
1248: /* Geometry */
1249: PetscFEGeom *cgeomFEM;
1250: DM dmGrad;
1251: PetscQuadrature affineQuad = NULL;
1252: Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1253: PetscFVCellGeom *cgeomFVM;
1254: const PetscScalar *lgrad;
1255: PetscInt maxDegree;
1256: DMField coordField;
1257: IS cellIS;
1258: PetscErrorCode ierr;
1261: DMGetDS(dm, &prob);
1262: DMGetDimension(dm, &dim);
1263: DMGetSection(dm, §ion);
1264: PetscSectionGetNumFields(section, &Nf);
1265: /* Determine which discretizations we have */
1266: for (f = 0; f < Nf; ++f) {
1267: PetscObject obj;
1268: PetscClassId id;
1270: PetscDSGetDiscretization(prob, f, &obj);
1271: PetscObjectGetClassId(obj, &id);
1272: if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1273: }
1274: /* Get local solution with boundary values */
1275: DMGetLocalVector(dm, &locX);
1276: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1277: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1278: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1279: /* Read DS information */
1280: PetscDSGetTotalDimension(prob, &totDim);
1281: PetscDSGetComponentOffsets(prob, &uOff);
1282: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1283: ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1284: PetscDSGetConstants(prob, &numConstants, &constants);
1285: /* Read Auxiliary DS information */
1286: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1287: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1288: if (dmAux) {
1289: DMGetDS(dmAux, &probAux);
1290: PetscDSGetNumFields(probAux, &NfAux);
1291: DMGetSection(dmAux, §ionAux);
1292: PetscDSGetTotalDimension(probAux, &totDimAux);
1293: PetscDSGetComponentOffsets(probAux, &aOff);
1294: }
1295: /* Allocate data arrays */
1296: PetscCalloc1(numCells*totDim, &u);
1297: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1298: /* Read out geometry */
1299: DMGetCoordinateField(dm,&coordField);
1300: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1301: if (maxDegree <= 1) {
1302: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1303: if (affineQuad) {
1304: DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1305: }
1306: }
1307: if (useFVM) {
1308: PetscFV fv = NULL;
1309: Vec grad;
1310: PetscInt fStart, fEnd;
1311: PetscBool compGrad;
1313: for (f = 0; f < Nf; ++f) {
1314: PetscObject obj;
1315: PetscClassId id;
1317: PetscDSGetDiscretization(prob, f, &obj);
1318: PetscObjectGetClassId(obj, &id);
1319: if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1320: }
1321: PetscFVGetComputeGradients(fv, &compGrad);
1322: PetscFVSetComputeGradients(fv, PETSC_TRUE);
1323: DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1324: DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1325: PetscFVSetComputeGradients(fv, compGrad);
1326: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1327: /* Reconstruct and limit cell gradients */
1328: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1329: DMGetGlobalVector(dmGrad, &grad);
1330: DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1331: /* Communicate gradient values */
1332: DMGetLocalVector(dmGrad, &locGrad);
1333: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1334: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1335: DMRestoreGlobalVector(dmGrad, &grad);
1336: /* Handle non-essential (e.g. outflow) boundary values */
1337: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1338: VecGetArrayRead(locGrad, &lgrad);
1339: }
1340: /* Read out data from inputs */
1341: for (c = cStart; c < cEnd; ++c) {
1342: PetscScalar *x = NULL;
1343: PetscInt i;
1345: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1346: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1347: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1348: if (dmAux) {
1349: DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1350: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1351: DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1352: }
1353: }
1354: /* Do integration for each field */
1355: for (f = 0; f < Nf; ++f) {
1356: PetscObject obj;
1357: PetscClassId id;
1358: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1360: PetscDSGetDiscretization(prob, f, &obj);
1361: PetscObjectGetClassId(obj, &id);
1362: if (id == PETSCFE_CLASSID) {
1363: PetscFE fe = (PetscFE) obj;
1364: PetscQuadrature q;
1365: PetscFEGeom *chunkGeom = NULL;
1366: PetscInt Nq, Nb;
1368: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1369: PetscFEGetQuadrature(fe, &q);
1370: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1371: PetscFEGetDimension(fe, &Nb);
1372: blockSize = Nb*Nq;
1373: batchSize = numBlocks * blockSize;
1374: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1375: numChunks = numCells / (numBatches*batchSize);
1376: Ne = numChunks*numBatches*batchSize;
1377: Nr = numCells % (numBatches*batchSize);
1378: offset = numCells - Nr;
1379: if (!affineQuad) {
1380: DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
1381: }
1382: PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1383: PetscFEIntegrate(fe, prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
1384: PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
1385: PetscFEIntegrate(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1386: PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
1387: if (!affineQuad) {
1388: PetscFEGeomDestroy(&cgeomFEM);
1389: }
1390: } else if (id == PETSCFV_CLASSID) {
1391: PetscInt foff;
1392: PetscPointFunc obj_func;
1393: PetscScalar lint;
1395: PetscDSGetObjective(prob, f, &obj_func);
1396: PetscDSGetFieldOffset(prob, f, &foff);
1397: if (obj_func) {
1398: for (c = 0; c < numCells; ++c) {
1399: PetscScalar *u_x;
1401: DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1402: 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);
1403: cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1404: }
1405: }
1406: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1407: }
1408: /* Cleanup data arrays */
1409: if (useFVM) {
1410: VecRestoreArrayRead(locGrad, &lgrad);
1411: VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1412: DMRestoreLocalVector(dmGrad, &locGrad);
1413: VecDestroy(&faceGeometryFVM);
1414: VecDestroy(&cellGeometryFVM);
1415: DMDestroy(&dmGrad);
1416: }
1417: if (dmAux) {PetscFree(a);}
1418: PetscFree(u);
1419: /* Cleanup */
1420: if (affineQuad) {
1421: PetscFEGeomDestroy(&cgeomFEM);
1422: }
1423: PetscQuadratureDestroy(&affineQuad);
1424: ISDestroy(&cellIS);
1425: DMRestoreLocalVector(dm, &locX);
1426: return(0);
1427: }
1429: /*@
1430: DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user
1432: Input Parameters:
1433: + dm - The mesh
1434: . X - Global input vector
1435: - user - The user context
1437: Output Parameter:
1438: . integral - Integral for each field
1440: Level: developer
1442: .seealso: DMPlexComputeResidualFEM()
1443: @*/
1444: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1445: {
1446: DM_Plex *mesh = (DM_Plex *) dm->data;
1447: PetscScalar *cintegral, *lintegral;
1448: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1455: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1456: DMGetNumFields(dm, &Nf);
1457: DMPlexGetVTKCellHeight(dm, &cellHeight);
1458: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1459: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1460: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1461: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1462: PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1463: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1464: /* Sum up values */
1465: for (cell = cStart; cell < cEnd; ++cell) {
1466: const PetscInt c = cell - cStart;
1468: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1469: for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1470: }
1471: MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1472: if (mesh->printFEM) {
1473: PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1474: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1475: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1476: }
1477: PetscFree2(lintegral, cintegral);
1478: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1479: return(0);
1480: }
1482: /*@
1483: DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user
1485: Input Parameters:
1486: + dm - The mesh
1487: . X - Global input vector
1488: - user - The user context
1490: Output Parameter:
1491: . integral - Cellwise integrals for each field
1493: Level: developer
1495: .seealso: DMPlexComputeResidualFEM()
1496: @*/
1497: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
1498: {
1499: DM_Plex *mesh = (DM_Plex *) dm->data;
1500: DM dmF;
1501: PetscSection sectionF;
1502: PetscScalar *cintegral, *af;
1503: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1510: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1511: DMGetNumFields(dm, &Nf);
1512: DMPlexGetVTKCellHeight(dm, &cellHeight);
1513: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1514: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1515: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1516: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1517: PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
1518: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1519: /* Put values in F*/
1520: VecGetDM(F, &dmF);
1521: DMGetSection(dmF, §ionF);
1522: VecGetArray(F, &af);
1523: for (cell = cStart; cell < cEnd; ++cell) {
1524: const PetscInt c = cell - cStart;
1525: PetscInt dof, off;
1527: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1528: PetscSectionGetDof(sectionF, cell, &dof);
1529: PetscSectionGetOffset(sectionF, cell, &off);
1530: if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
1531: for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
1532: }
1533: VecRestoreArray(F, &af);
1534: PetscFree(cintegral);
1535: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1536: return(0);
1537: }
1539: static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS,
1540: void (*func)(PetscInt, PetscInt, PetscInt,
1541: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1542: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1543: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1544: PetscScalar *fintegral, void *user)
1545: {
1546: DM plex = NULL, plexA = NULL;
1547: PetscDS prob, probAux = NULL;
1548: PetscSection section, sectionAux = NULL;
1549: Vec locA = NULL;
1550: DMField coordField;
1551: PetscInt Nf, totDim, *uOff, *uOff_x;
1552: PetscInt NfAux = 0, totDimAux = 0, *aOff = NULL;
1553: PetscScalar *u, *a = NULL;
1554: const PetscScalar *constants;
1555: PetscInt numConstants, f;
1556: PetscErrorCode ierr;
1559: DMGetCoordinateField(dm, &coordField);
1560: DMConvert(dm, DMPLEX, &plex);
1561: DMGetDS(dm, &prob);
1562: DMGetDefaultSection(dm, §ion);
1563: PetscSectionGetNumFields(section, &Nf);
1564: /* Determine which discretizations we have */
1565: for (f = 0; f < Nf; ++f) {
1566: PetscObject obj;
1567: PetscClassId id;
1569: PetscDSGetDiscretization(prob, f, &obj);
1570: PetscObjectGetClassId(obj, &id);
1571: if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f);
1572: }
1573: /* Read DS information */
1574: PetscDSGetTotalDimension(prob, &totDim);
1575: PetscDSGetComponentOffsets(prob, &uOff);
1576: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1577: PetscDSGetConstants(prob, &numConstants, &constants);
1578: /* Read Auxiliary DS information */
1579: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1580: if (locA) {
1581: DM dmAux;
1583: VecGetDM(locA, &dmAux);
1584: DMConvert(dmAux, DMPLEX, &plexA);
1585: DMGetDS(dmAux, &probAux);
1586: PetscDSGetNumFields(probAux, &NfAux);
1587: DMGetDefaultSection(dmAux, §ionAux);
1588: PetscDSGetTotalDimension(probAux, &totDimAux);
1589: PetscDSGetComponentOffsets(probAux, &aOff);
1590: }
1591: /* Integrate over points */
1592: {
1593: PetscFEGeom *fgeom, *chunkGeom = NULL;
1594: PetscInt maxDegree;
1595: PetscQuadrature qGeom = NULL;
1596: const PetscInt *points;
1597: PetscInt numFaces, face, Nq, field;
1598: PetscInt numChunks, chunkSize, chunk, Nr, offset;
1600: ISGetLocalSize(pointIS, &numFaces);
1601: ISGetIndices(pointIS, &points);
1602: PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);
1603: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
1604: for (field = 0; field < Nf; ++field) {
1605: PetscFE fe;
1607: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1608: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);}
1609: if (!qGeom) {
1610: PetscFEGetFaceQuadrature(fe, &qGeom);
1611: PetscObjectReference((PetscObject) qGeom);
1612: }
1613: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
1614: DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
1615: for (face = 0; face < numFaces; ++face) {
1616: const PetscInt point = points[face], *support, *cone;
1617: PetscScalar *x = NULL;
1618: PetscInt i, coneSize, faceLoc;
1620: DMPlexGetSupport(dm, point, &support);
1621: DMPlexGetConeSize(dm, support[0], &coneSize);
1622: DMPlexGetCone(dm, support[0], &cone);
1623: for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
1624: if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]);
1625: fgeom->face[face][0] = faceLoc;
1626: DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
1627: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
1628: DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
1629: if (locA) {
1630: PetscInt subp;
1631: DMPlexGetSubpoint(plexA, support[0], &subp);
1632: DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
1633: for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i];
1634: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
1635: }
1636: }
1637: /* Get blocking */
1638: {
1639: PetscQuadrature q;
1640: PetscInt numBatches, batchSize, numBlocks, blockSize;
1641: PetscInt Nq, Nb;
1643: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1644: PetscFEGetQuadrature(fe, &q);
1645: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1646: PetscFEGetDimension(fe, &Nb);
1647: blockSize = Nb*Nq;
1648: batchSize = numBlocks * blockSize;
1649: chunkSize = numBatches*batchSize;
1650: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1651: numChunks = numFaces / chunkSize;
1652: Nr = numFaces % chunkSize;
1653: offset = numFaces - Nr;
1654: }
1655: /* Do integration for each field */
1656: for (chunk = 0; chunk < numChunks; ++chunk) {
1657: PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);
1658: PetscFEIntegrateBd(fe, prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);
1659: PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
1660: }
1661: PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);
1662: PetscFEIntegrateBd(fe, prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);
1663: PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);
1664: /* Cleanup data arrays */
1665: DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
1666: PetscQuadratureDestroy(&qGeom);
1667: PetscFree2(u, a);
1668: ISRestoreIndices(pointIS, &points);
1669: }
1670: }
1671: if (plex) {DMDestroy(&plex);}
1672: if (plexA) {DMDestroy(&plexA);}
1673: return(0);
1674: }
1676: /*@
1677: DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user
1679: Input Parameters:
1680: + dm - The mesh
1681: . X - Global input vector
1682: . label - The boundary DMLabel
1683: . numVals - The number of label values to use, or PETSC_DETERMINE for all values
1684: . vals - The label values to use, or PETSC_NULL for all values
1685: . func = The function to integrate along the boundary
1686: - user - The user context
1688: Output Parameter:
1689: . integral - Integral for each field
1691: Level: developer
1693: .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM()
1694: @*/
1695: PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[],
1696: void (*func)(PetscInt, PetscInt, PetscInt,
1697: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1698: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1699: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1700: PetscScalar *integral, void *user)
1701: {
1702: Vec locX;
1703: PetscSection section;
1704: DMLabel depthLabel;
1705: IS facetIS;
1706: PetscInt dim, Nf, f, v;
1715: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1716: DMPlexGetDepthLabel(dm, &depthLabel);
1717: DMGetDimension(dm, &dim);
1718: DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
1719: DMGetDefaultSection(dm, §ion);
1720: PetscSectionGetNumFields(section, &Nf);
1721: /* Get local solution with boundary values */
1722: DMGetLocalVector(dm, &locX);
1723: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1724: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1725: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1726: /* Loop over label values */
1727: PetscMemzero(integral, Nf * sizeof(PetscScalar));
1728: for (v = 0; v < numVals; ++v) {
1729: IS pointIS;
1730: PetscInt numFaces, face;
1731: PetscScalar *fintegral;
1733: DMLabelGetStratumIS(label, vals[v], &pointIS);
1734: if (!pointIS) continue; /* No points with that id on this process */
1735: {
1736: IS isectIS;
1738: /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
1739: ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);
1740: ISDestroy(&pointIS);
1741: pointIS = isectIS;
1742: }
1743: ISGetLocalSize(pointIS, &numFaces);
1744: PetscCalloc1(numFaces*Nf, &fintegral);
1745: DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);
1746: /* Sum point contributions into integral */
1747: for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f];
1748: PetscFree(fintegral);
1749: ISDestroy(&pointIS);
1750: }
1751: DMRestoreLocalVector(dm, &locX);
1752: ISDestroy(&facetIS);
1753: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1754: return(0);
1755: }
1757: /*@
1758: DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1760: Input Parameters:
1761: + dmf - The fine mesh
1762: . dmc - The coarse mesh
1763: - user - The user context
1765: Output Parameter:
1766: . In - The interpolation matrix
1768: Level: developer
1770: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1771: @*/
1772: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1773: {
1774: DM_Plex *mesh = (DM_Plex *) dmc->data;
1775: const char *name = "Interpolator";
1776: PetscDS prob;
1777: PetscFE *feRef;
1778: PetscFV *fvRef;
1779: PetscSection fsection, fglobalSection;
1780: PetscSection csection, cglobalSection;
1781: PetscScalar *elemMat;
1782: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1783: PetscInt cTotDim, rTotDim = 0;
1784: PetscErrorCode ierr;
1787: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1788: DMGetDimension(dmf, &dim);
1789: DMGetSection(dmf, &fsection);
1790: DMGetGlobalSection(dmf, &fglobalSection);
1791: DMGetSection(dmc, &csection);
1792: DMGetGlobalSection(dmc, &cglobalSection);
1793: PetscSectionGetNumFields(fsection, &Nf);
1794: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1795: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1796: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1797: DMGetDS(dmf, &prob);
1798: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1799: for (f = 0; f < Nf; ++f) {
1800: PetscObject obj;
1801: PetscClassId id;
1802: PetscInt rNb = 0, Nc = 0;
1804: PetscDSGetDiscretization(prob, f, &obj);
1805: PetscObjectGetClassId(obj, &id);
1806: if (id == PETSCFE_CLASSID) {
1807: PetscFE fe = (PetscFE) obj;
1809: PetscFERefine(fe, &feRef[f]);
1810: PetscFEGetDimension(feRef[f], &rNb);
1811: PetscFEGetNumComponents(fe, &Nc);
1812: } else if (id == PETSCFV_CLASSID) {
1813: PetscFV fv = (PetscFV) obj;
1814: PetscDualSpace Q;
1816: PetscFVRefine(fv, &fvRef[f]);
1817: PetscFVGetDualSpace(fvRef[f], &Q);
1818: PetscDualSpaceGetDimension(Q, &rNb);
1819: PetscFVGetNumComponents(fv, &Nc);
1820: }
1821: rTotDim += rNb;
1822: }
1823: PetscDSGetTotalDimension(prob, &cTotDim);
1824: PetscMalloc1(rTotDim*cTotDim,&elemMat);
1825: PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1826: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1827: PetscDualSpace Qref;
1828: PetscQuadrature f;
1829: const PetscReal *qpoints, *qweights;
1830: PetscReal *points;
1831: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
1833: /* Compose points from all dual basis functionals */
1834: if (feRef[fieldI]) {
1835: PetscFEGetDualSpace(feRef[fieldI], &Qref);
1836: PetscFEGetNumComponents(feRef[fieldI], &Nc);
1837: } else {
1838: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1839: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1840: }
1841: PetscDualSpaceGetDimension(Qref, &fpdim);
1842: for (i = 0; i < fpdim; ++i) {
1843: PetscDualSpaceGetFunctional(Qref, i, &f);
1844: PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1845: npoints += Np;
1846: }
1847: PetscMalloc1(npoints*dim,&points);
1848: for (i = 0, k = 0; i < fpdim; ++i) {
1849: PetscDualSpaceGetFunctional(Qref, i, &f);
1850: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1851: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1852: }
1854: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1855: PetscObject obj;
1856: PetscClassId id;
1857: PetscReal *B;
1858: PetscInt NcJ = 0, cpdim = 0, j, qNc;
1860: PetscDSGetDiscretization(prob, fieldJ, &obj);
1861: PetscObjectGetClassId(obj, &id);
1862: if (id == PETSCFE_CLASSID) {
1863: PetscFE fe = (PetscFE) obj;
1865: /* Evaluate basis at points */
1866: PetscFEGetNumComponents(fe, &NcJ);
1867: PetscFEGetDimension(fe, &cpdim);
1868: /* For now, fields only interpolate themselves */
1869: if (fieldI == fieldJ) {
1870: 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);
1871: PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1872: for (i = 0, k = 0; i < fpdim; ++i) {
1873: PetscDualSpaceGetFunctional(Qref, i, &f);
1874: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1875: 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);
1876: for (p = 0; p < Np; ++p, ++k) {
1877: for (j = 0; j < cpdim; ++j) {
1878: /*
1879: cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
1880: offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields
1881: fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
1882: qNC, Nc, Ncj, c: Number of components in this field
1883: Np, p: Number of quad points in the fine grid functional i
1884: k: i*Np + p, overall point number for the interpolation
1885: */
1886: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1887: }
1888: }
1889: }
1890: PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1891: }
1892: } else if (id == PETSCFV_CLASSID) {
1893: PetscFV fv = (PetscFV) obj;
1895: /* Evaluate constant function at points */
1896: PetscFVGetNumComponents(fv, &NcJ);
1897: cpdim = 1;
1898: /* For now, fields only interpolate themselves */
1899: if (fieldI == fieldJ) {
1900: 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);
1901: for (i = 0, k = 0; i < fpdim; ++i) {
1902: PetscDualSpaceGetFunctional(Qref, i, &f);
1903: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1904: 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);
1905: for (p = 0; p < Np; ++p, ++k) {
1906: for (j = 0; j < cpdim; ++j) {
1907: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c];
1908: }
1909: }
1910: }
1911: }
1912: }
1913: offsetJ += cpdim;
1914: }
1915: offsetI += fpdim;
1916: PetscFree(points);
1917: }
1918: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1919: /* Preallocate matrix */
1920: {
1921: Mat preallocator;
1922: PetscScalar *vals;
1923: PetscInt *cellCIndices, *cellFIndices;
1924: PetscInt locRows, locCols, cell;
1926: MatGetLocalSize(In, &locRows, &locCols);
1927: MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1928: MatSetType(preallocator, MATPREALLOCATOR);
1929: MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1930: MatSetUp(preallocator);
1931: PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1932: for (cell = cStart; cell < cEnd; ++cell) {
1933: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1934: MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1935: }
1936: PetscFree3(vals,cellCIndices,cellFIndices);
1937: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1938: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1939: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1940: MatDestroy(&preallocator);
1941: }
1942: /* Fill matrix */
1943: MatZeroEntries(In);
1944: for (c = cStart; c < cEnd; ++c) {
1945: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1946: }
1947: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1948: PetscFree2(feRef,fvRef);
1949: PetscFree(elemMat);
1950: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1951: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1952: if (mesh->printFEM) {
1953: PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);
1954: MatChop(In, 1.0e-10);
1955: MatView(In, NULL);
1956: }
1957: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1958: return(0);
1959: }
1961: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
1962: {
1963: SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
1964: }
1966: /*@
1967: DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.
1969: Input Parameters:
1970: + dmf - The fine mesh
1971: . dmc - The coarse mesh
1972: - user - The user context
1974: Output Parameter:
1975: . In - The interpolation matrix
1977: Level: developer
1979: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1980: @*/
1981: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1982: {
1983: DM_Plex *mesh = (DM_Plex *) dmf->data;
1984: const char *name = "Interpolator";
1985: PetscDS prob;
1986: PetscSection fsection, csection, globalFSection, globalCSection;
1987: PetscHSetIJ ht;
1988: PetscLayout rLayout;
1989: PetscInt *dnz, *onz;
1990: PetscInt locRows, rStart, rEnd;
1991: PetscReal *x, *v0, *J, *invJ, detJ;
1992: PetscReal *v0c, *Jc, *invJc, detJc;
1993: PetscScalar *elemMat;
1994: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
1998: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1999: DMGetCoordinateDim(dmc, &dim);
2000: DMGetDS(dmc, &prob);
2001: PetscDSGetRefCoordArrays(prob, &x, NULL);
2002: PetscDSGetNumFields(prob, &Nf);
2003: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2004: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2005: DMGetSection(dmf, &fsection);
2006: DMGetGlobalSection(dmf, &globalFSection);
2007: DMGetSection(dmc, &csection);
2008: DMGetGlobalSection(dmc, &globalCSection);
2009: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2010: PetscDSGetTotalDimension(prob, &totDim);
2011: PetscMalloc1(totDim, &elemMat);
2013: MatGetLocalSize(In, &locRows, NULL);
2014: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
2015: PetscLayoutSetLocalSize(rLayout, locRows);
2016: PetscLayoutSetBlockSize(rLayout, 1);
2017: PetscLayoutSetUp(rLayout);
2018: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2019: PetscLayoutDestroy(&rLayout);
2020: PetscCalloc2(locRows,&dnz,locRows,&onz);
2021: PetscHSetIJCreate(&ht);
2022: for (field = 0; field < Nf; ++field) {
2023: PetscObject obj;
2024: PetscClassId id;
2025: PetscDualSpace Q = NULL;
2026: PetscQuadrature f;
2027: const PetscReal *qpoints;
2028: PetscInt Nc, Np, fpdim, i, d;
2030: PetscDSGetDiscretization(prob, field, &obj);
2031: PetscObjectGetClassId(obj, &id);
2032: if (id == PETSCFE_CLASSID) {
2033: PetscFE fe = (PetscFE) obj;
2035: PetscFEGetDualSpace(fe, &Q);
2036: PetscFEGetNumComponents(fe, &Nc);
2037: } else if (id == PETSCFV_CLASSID) {
2038: PetscFV fv = (PetscFV) obj;
2040: PetscFVGetDualSpace(fv, &Q);
2041: Nc = 1;
2042: }
2043: PetscDualSpaceGetDimension(Q, &fpdim);
2044: /* For each fine grid cell */
2045: for (cell = cStart; cell < cEnd; ++cell) {
2046: PetscInt *findices, *cindices;
2047: PetscInt numFIndices, numCIndices;
2049: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2050: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2051: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2052: for (i = 0; i < fpdim; ++i) {
2053: Vec pointVec;
2054: PetscScalar *pV;
2055: PetscSF coarseCellSF = NULL;
2056: const PetscSFNode *coarseCells;
2057: PetscInt numCoarseCells, q, c;
2059: /* Get points from the dual basis functional quadrature */
2060: PetscDualSpaceGetFunctional(Q, i, &f);
2061: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2062: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2063: VecSetBlockSize(pointVec, dim);
2064: VecGetArray(pointVec, &pV);
2065: for (q = 0; q < Np; ++q) {
2066: const PetscReal xi0[3] = {-1., -1., -1.};
2068: /* Transform point to real space */
2069: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2070: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2071: }
2072: VecRestoreArray(pointVec, &pV);
2073: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2074: /* OPT: Pack all quad points from fine cell */
2075: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2076: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2077: /* Update preallocation info */
2078: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2079: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2080: {
2081: PetscHashIJKey key;
2082: PetscBool missing;
2084: key.i = findices[i];
2085: if (key.i >= 0) {
2086: /* Get indices for coarse elements */
2087: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2088: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2089: for (c = 0; c < numCIndices; ++c) {
2090: key.j = cindices[c];
2091: if (key.j < 0) continue;
2092: PetscHSetIJQueryAdd(ht, key, &missing);
2093: if (missing) {
2094: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2095: else ++onz[key.i-rStart];
2096: }
2097: }
2098: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2099: }
2100: }
2101: }
2102: PetscSFDestroy(&coarseCellSF);
2103: VecDestroy(&pointVec);
2104: }
2105: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2106: }
2107: }
2108: PetscHSetIJDestroy(&ht);
2109: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
2110: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2111: PetscFree2(dnz,onz);
2112: for (field = 0; field < Nf; ++field) {
2113: PetscObject obj;
2114: PetscClassId id;
2115: PetscDualSpace Q = NULL;
2116: PetscQuadrature f;
2117: const PetscReal *qpoints, *qweights;
2118: PetscInt Nc, qNc, Np, fpdim, i, d;
2120: PetscDSGetDiscretization(prob, field, &obj);
2121: PetscObjectGetClassId(obj, &id);
2122: if (id == PETSCFE_CLASSID) {
2123: PetscFE fe = (PetscFE) obj;
2125: PetscFEGetDualSpace(fe, &Q);
2126: PetscFEGetNumComponents(fe, &Nc);
2127: } else if (id == PETSCFV_CLASSID) {
2128: PetscFV fv = (PetscFV) obj;
2130: PetscFVGetDualSpace(fv, &Q);
2131: Nc = 1;
2132: } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field);
2133: PetscDualSpaceGetDimension(Q, &fpdim);
2134: /* For each fine grid cell */
2135: for (cell = cStart; cell < cEnd; ++cell) {
2136: PetscInt *findices, *cindices;
2137: PetscInt numFIndices, numCIndices;
2139: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2140: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2141: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2142: for (i = 0; i < fpdim; ++i) {
2143: Vec pointVec;
2144: PetscScalar *pV;
2145: PetscSF coarseCellSF = NULL;
2146: const PetscSFNode *coarseCells;
2147: PetscInt numCoarseCells, cpdim, q, c, j;
2149: /* Get points from the dual basis functional quadrature */
2150: PetscDualSpaceGetFunctional(Q, i, &f);
2151: PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
2152: 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);
2153: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2154: VecSetBlockSize(pointVec, dim);
2155: VecGetArray(pointVec, &pV);
2156: for (q = 0; q < Np; ++q) {
2157: const PetscReal xi0[3] = {-1., -1., -1.};
2159: /* Transform point to real space */
2160: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2161: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2162: }
2163: VecRestoreArray(pointVec, &pV);
2164: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2165: /* OPT: Read this out from preallocation information */
2166: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2167: /* Update preallocation info */
2168: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2169: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2170: VecGetArray(pointVec, &pV);
2171: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2172: PetscReal pVReal[3];
2173: const PetscReal xi0[3] = {-1., -1., -1.};
2175: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2176: /* Transform points from real space to coarse reference space */
2177: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2178: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2179: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
2181: if (id == PETSCFE_CLASSID) {
2182: PetscFE fe = (PetscFE) obj;
2183: PetscReal *B;
2185: /* Evaluate coarse basis on contained point */
2186: PetscFEGetDimension(fe, &cpdim);
2187: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2188: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2189: /* Get elemMat entries by multiplying by weight */
2190: for (j = 0; j < cpdim; ++j) {
2191: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
2192: }
2193: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2194: } else {
2195: cpdim = 1;
2196: for (j = 0; j < cpdim; ++j) {
2197: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
2198: }
2199: }
2200: /* Update interpolator */
2201: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2202: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2203: MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
2204: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2205: }
2206: VecRestoreArray(pointVec, &pV);
2207: PetscSFDestroy(&coarseCellSF);
2208: VecDestroy(&pointVec);
2209: }
2210: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2211: }
2212: }
2213: PetscFree3(v0,J,invJ);
2214: PetscFree3(v0c,Jc,invJc);
2215: PetscFree(elemMat);
2216: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2217: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2218: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2219: return(0);
2220: }
2222: /*@
2223: DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM.
2225: Input Parameters:
2226: + dmf - The fine mesh
2227: . dmc - The coarse mesh
2228: - user - The user context
2230: Output Parameter:
2231: . mass - The mass matrix
2233: Level: developer
2235: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2236: @*/
2237: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
2238: {
2239: DM_Plex *mesh = (DM_Plex *) dmf->data;
2240: const char *name = "Mass Matrix";
2241: PetscDS prob;
2242: PetscSection fsection, csection, globalFSection, globalCSection;
2243: PetscHSetIJ ht;
2244: PetscLayout rLayout;
2245: PetscInt *dnz, *onz;
2246: PetscInt locRows, rStart, rEnd;
2247: PetscReal *x, *v0, *J, *invJ, detJ;
2248: PetscReal *v0c, *Jc, *invJc, detJc;
2249: PetscScalar *elemMat;
2250: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
2254: DMGetCoordinateDim(dmc, &dim);
2255: DMGetDS(dmc, &prob);
2256: PetscDSGetRefCoordArrays(prob, &x, NULL);
2257: PetscDSGetNumFields(prob, &Nf);
2258: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2259: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2260: DMGetSection(dmf, &fsection);
2261: DMGetGlobalSection(dmf, &globalFSection);
2262: DMGetSection(dmc, &csection);
2263: DMGetGlobalSection(dmc, &globalCSection);
2264: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2265: PetscDSGetTotalDimension(prob, &totDim);
2266: PetscMalloc1(totDim, &elemMat);
2268: MatGetLocalSize(mass, &locRows, NULL);
2269: PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
2270: PetscLayoutSetLocalSize(rLayout, locRows);
2271: PetscLayoutSetBlockSize(rLayout, 1);
2272: PetscLayoutSetUp(rLayout);
2273: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2274: PetscLayoutDestroy(&rLayout);
2275: PetscCalloc2(locRows,&dnz,locRows,&onz);
2276: PetscHSetIJCreate(&ht);
2277: for (field = 0; field < Nf; ++field) {
2278: PetscObject obj;
2279: PetscClassId id;
2280: PetscQuadrature quad;
2281: const PetscReal *qpoints;
2282: PetscInt Nq, Nc, i, d;
2284: PetscDSGetDiscretization(prob, field, &obj);
2285: PetscObjectGetClassId(obj, &id);
2286: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
2287: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2288: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
2289: /* For each fine grid cell */
2290: for (cell = cStart; cell < cEnd; ++cell) {
2291: Vec pointVec;
2292: PetscScalar *pV;
2293: PetscSF coarseCellSF = NULL;
2294: const PetscSFNode *coarseCells;
2295: PetscInt numCoarseCells, q, c;
2296: PetscInt *findices, *cindices;
2297: PetscInt numFIndices, numCIndices;
2299: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2300: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2301: /* Get points from the quadrature */
2302: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2303: VecSetBlockSize(pointVec, dim);
2304: VecGetArray(pointVec, &pV);
2305: for (q = 0; q < Nq; ++q) {
2306: const PetscReal xi0[3] = {-1., -1., -1.};
2308: /* Transform point to real space */
2309: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2310: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2311: }
2312: VecRestoreArray(pointVec, &pV);
2313: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2314: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2315: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2316: /* Update preallocation info */
2317: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2318: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2319: {
2320: PetscHashIJKey key;
2321: PetscBool missing;
2323: for (i = 0; i < numFIndices; ++i) {
2324: key.i = findices[i];
2325: if (key.i >= 0) {
2326: /* Get indices for coarse elements */
2327: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2328: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2329: for (c = 0; c < numCIndices; ++c) {
2330: key.j = cindices[c];
2331: if (key.j < 0) continue;
2332: PetscHSetIJQueryAdd(ht, key, &missing);
2333: if (missing) {
2334: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2335: else ++onz[key.i-rStart];
2336: }
2337: }
2338: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2339: }
2340: }
2341: }
2342: }
2343: PetscSFDestroy(&coarseCellSF);
2344: VecDestroy(&pointVec);
2345: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2346: }
2347: }
2348: PetscHSetIJDestroy(&ht);
2349: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2350: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2351: PetscFree2(dnz,onz);
2352: for (field = 0; field < Nf; ++field) {
2353: PetscObject obj;
2354: PetscClassId id;
2355: PetscQuadrature quad;
2356: PetscReal *Bfine;
2357: const PetscReal *qpoints, *qweights;
2358: PetscInt Nq, Nc, i, d;
2360: PetscDSGetDiscretization(prob, field, &obj);
2361: PetscObjectGetClassId(obj, &id);
2362: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);}
2363: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2364: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
2365: /* For each fine grid cell */
2366: for (cell = cStart; cell < cEnd; ++cell) {
2367: Vec pointVec;
2368: PetscScalar *pV;
2369: PetscSF coarseCellSF = NULL;
2370: const PetscSFNode *coarseCells;
2371: PetscInt numCoarseCells, cpdim, q, c, j;
2372: PetscInt *findices, *cindices;
2373: PetscInt numFIndices, numCIndices;
2375: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2376: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2377: /* Get points from the quadrature */
2378: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2379: VecSetBlockSize(pointVec, dim);
2380: VecGetArray(pointVec, &pV);
2381: for (q = 0; q < Nq; ++q) {
2382: const PetscReal xi0[3] = {-1., -1., -1.};
2384: /* Transform point to real space */
2385: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2386: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2387: }
2388: VecRestoreArray(pointVec, &pV);
2389: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2390: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2391: /* Update matrix */
2392: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2393: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2394: VecGetArray(pointVec, &pV);
2395: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2396: PetscReal pVReal[3];
2397: const PetscReal xi0[3] = {-1., -1., -1.};
2400: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2401: /* Transform points from real space to coarse reference space */
2402: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2403: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2404: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
2406: if (id == PETSCFE_CLASSID) {
2407: PetscFE fe = (PetscFE) obj;
2408: PetscReal *B;
2410: /* Evaluate coarse basis on contained point */
2411: PetscFEGetDimension(fe, &cpdim);
2412: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2413: /* Get elemMat entries by multiplying by weight */
2414: for (i = 0; i < numFIndices; ++i) {
2415: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2416: for (j = 0; j < cpdim; ++j) {
2417: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2418: }
2419: /* Update interpolator */
2420: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2421: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2422: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2423: }
2424: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2425: } else {
2426: cpdim = 1;
2427: for (i = 0; i < numFIndices; ++i) {
2428: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2429: for (j = 0; j < cpdim; ++j) {
2430: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2431: }
2432: /* Update interpolator */
2433: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2434: PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);
2435: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2436: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2437: }
2438: }
2439: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2440: }
2441: VecRestoreArray(pointVec, &pV);
2442: PetscSFDestroy(&coarseCellSF);
2443: VecDestroy(&pointVec);
2444: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2445: }
2446: }
2447: PetscFree3(v0,J,invJ);
2448: PetscFree3(v0c,Jc,invJc);
2449: PetscFree(elemMat);
2450: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2451: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2452: return(0);
2453: }
2455: /*@
2456: DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns
2458: Input Parameters:
2459: + dmc - The coarse mesh
2460: - dmf - The fine mesh
2461: - user - The user context
2463: Output Parameter:
2464: . sc - The mapping
2466: Level: developer
2468: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2469: @*/
2470: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
2471: {
2472: PetscDS prob;
2473: PetscFE *feRef;
2474: PetscFV *fvRef;
2475: Vec fv, cv;
2476: IS fis, cis;
2477: PetscSection fsection, fglobalSection, csection, cglobalSection;
2478: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
2479: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
2480: PetscBool *needAvg;
2484: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2485: DMGetDimension(dmf, &dim);
2486: DMGetSection(dmf, &fsection);
2487: DMGetGlobalSection(dmf, &fglobalSection);
2488: DMGetSection(dmc, &csection);
2489: DMGetGlobalSection(dmc, &cglobalSection);
2490: PetscSectionGetNumFields(fsection, &Nf);
2491: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2492: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2493: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2494: DMGetDS(dmc, &prob);
2495: PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);
2496: for (f = 0; f < Nf; ++f) {
2497: PetscObject obj;
2498: PetscClassId id;
2499: PetscInt fNb = 0, Nc = 0;
2501: PetscDSGetDiscretization(prob, f, &obj);
2502: PetscObjectGetClassId(obj, &id);
2503: if (id == PETSCFE_CLASSID) {
2504: PetscFE fe = (PetscFE) obj;
2505: PetscSpace sp;
2506: PetscInt maxDegree;
2508: PetscFERefine(fe, &feRef[f]);
2509: PetscFEGetDimension(feRef[f], &fNb);
2510: PetscFEGetNumComponents(fe, &Nc);
2511: PetscFEGetBasisSpace(fe, &sp);
2512: PetscSpaceGetDegree(sp, NULL, &maxDegree);
2513: if (!maxDegree) needAvg[f] = PETSC_TRUE;
2514: } else if (id == PETSCFV_CLASSID) {
2515: PetscFV fv = (PetscFV) obj;
2516: PetscDualSpace Q;
2518: PetscFVRefine(fv, &fvRef[f]);
2519: PetscFVGetDualSpace(fvRef[f], &Q);
2520: PetscDualSpaceGetDimension(Q, &fNb);
2521: PetscFVGetNumComponents(fv, &Nc);
2522: needAvg[f] = PETSC_TRUE;
2523: }
2524: fTotDim += fNb;
2525: }
2526: PetscDSGetTotalDimension(prob, &cTotDim);
2527: PetscMalloc1(cTotDim,&cmap);
2528: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
2529: PetscFE feC;
2530: PetscFV fvC;
2531: PetscDualSpace QF, QC;
2532: PetscInt order = -1, NcF, NcC, fpdim, cpdim;
2534: if (feRef[field]) {
2535: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
2536: PetscFEGetNumComponents(feC, &NcC);
2537: PetscFEGetNumComponents(feRef[field], &NcF);
2538: PetscFEGetDualSpace(feRef[field], &QF);
2539: PetscDualSpaceGetOrder(QF, &order);
2540: PetscDualSpaceGetDimension(QF, &fpdim);
2541: PetscFEGetDualSpace(feC, &QC);
2542: PetscDualSpaceGetDimension(QC, &cpdim);
2543: } else {
2544: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
2545: PetscFVGetNumComponents(fvC, &NcC);
2546: PetscFVGetNumComponents(fvRef[field], &NcF);
2547: PetscFVGetDualSpace(fvRef[field], &QF);
2548: PetscDualSpaceGetDimension(QF, &fpdim);
2549: PetscFVGetDualSpace(fvC, &QC);
2550: PetscDualSpaceGetDimension(QC, &cpdim);
2551: }
2552: 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);
2553: for (c = 0; c < cpdim; ++c) {
2554: PetscQuadrature cfunc;
2555: const PetscReal *cqpoints, *cqweights;
2556: PetscInt NqcC, NpC;
2557: PetscBool found = PETSC_FALSE;
2559: PetscDualSpaceGetFunctional(QC, c, &cfunc);
2560: PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
2561: if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC);
2562: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
2563: for (f = 0; f < fpdim; ++f) {
2564: PetscQuadrature ffunc;
2565: const PetscReal *fqpoints, *fqweights;
2566: PetscReal sum = 0.0;
2567: PetscInt NqcF, NpF;
2569: PetscDualSpaceGetFunctional(QF, f, &ffunc);
2570: PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
2571: if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF);
2572: if (NpC != NpF) continue;
2573: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
2574: if (sum > 1.0e-9) continue;
2575: for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
2576: if (sum < 1.0e-9) continue;
2577: cmap[offsetC+c] = offsetF+f;
2578: found = PETSC_TRUE;
2579: break;
2580: }
2581: if (!found) {
2582: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
2583: if (fvRef[field] || (feRef[field] && order == 0)) {
2584: cmap[offsetC+c] = offsetF+0;
2585: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
2586: }
2587: }
2588: offsetC += cpdim;
2589: offsetF += fpdim;
2590: }
2591: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
2592: PetscFree3(feRef,fvRef,needAvg);
2594: DMGetGlobalVector(dmf, &fv);
2595: DMGetGlobalVector(dmc, &cv);
2596: VecGetOwnershipRange(cv, &startC, &endC);
2597: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
2598: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
2599: PetscMalloc1(m,&cindices);
2600: PetscMalloc1(m,&findices);
2601: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
2602: for (c = cStart; c < cEnd; ++c) {
2603: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
2604: for (d = 0; d < cTotDim; ++d) {
2605: if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
2606: 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]]);
2607: cindices[cellCIndices[d]-startC] = cellCIndices[d];
2608: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
2609: }
2610: }
2611: PetscFree(cmap);
2612: PetscFree2(cellCIndices,cellFIndices);
2614: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
2615: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
2616: VecScatterCreate(cv, cis, fv, fis, sc);
2617: ISDestroy(&cis);
2618: ISDestroy(&fis);
2619: DMRestoreGlobalVector(dmf, &fv);
2620: DMRestoreGlobalVector(dmc, &cv);
2621: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2622: return(0);
2623: }
2625: /*@C
2626: DMPlexGetCellFields - Retrieve the field values values for a chunk of cells
2628: Input Parameters:
2629: + dm - The DM
2630: . cellIS - The cells to include
2631: . locX - A local vector with the solution fields
2632: . locX_t - A local vector with solution field time derivatives, or NULL
2633: - locA - A local vector with auxiliary fields, or NULL
2635: Output Parameters:
2636: + u - The field coefficients
2637: . u_t - The fields derivative coefficients
2638: - a - The auxiliary field coefficients
2640: Level: developer
2642: .seealso: DMPlexGetFaceFields()
2643: @*/
2644: PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
2645: {
2646: DM plex, plexA = NULL;
2647: PetscSection section, sectionAux;
2648: PetscDS prob;
2649: const PetscInt *cells;
2650: PetscInt cStart, cEnd, numCells, totDim, totDimAux, c;
2651: PetscErrorCode ierr;
2661: DMPlexConvertPlex(dm, &plex, PETSC_FALSE);
2662: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
2663: DMGetSection(dm, §ion);
2664: DMGetCellDS(dm, cStart, &prob);
2665: PetscDSGetTotalDimension(prob, &totDim);
2666: if (locA) {
2667: DM dmAux;
2668: PetscDS probAux;
2670: VecGetDM(locA, &dmAux);
2671: DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);
2672: DMGetSection(dmAux, §ionAux);
2673: DMGetDS(dmAux, &probAux);
2674: PetscDSGetTotalDimension(probAux, &totDimAux);
2675: }
2676: numCells = cEnd - cStart;
2677: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);
2678: if (locX_t) {DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);} else {*u_t = NULL;}
2679: if (locA) {DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);} else {*a = NULL;}
2680: for (c = cStart; c < cEnd; ++c) {
2681: const PetscInt cell = cells ? cells[c] : c;
2682: const PetscInt cind = c - cStart;
2683: PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
2684: PetscInt i;
2686: DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);
2687: for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i];
2688: DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);
2689: if (locX_t) {
2690: DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);
2691: for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i];
2692: DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);
2693: }
2694: if (locA) {
2695: PetscInt subcell;
2696: DMPlexGetAuxiliaryPoint(plex, plexA, cell, &subcell);
2697: DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);
2698: for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i];
2699: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);
2700: }
2701: }
2702: DMDestroy(&plex);
2703: if (locA) {DMDestroy(&plexA);}
2704: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
2705: return(0);
2706: }
2708: /*@C
2709: DMPlexRestoreCellFields - Restore the field values values for a chunk of cells
2711: Input Parameters:
2712: + dm - The DM
2713: . cellIS - The cells to include
2714: . locX - A local vector with the solution fields
2715: . locX_t - A local vector with solution field time derivatives, or NULL
2716: - locA - A local vector with auxiliary fields, or NULL
2718: Output Parameters:
2719: + u - The field coefficients
2720: . u_t - The fields derivative coefficients
2721: - a - The auxiliary field coefficients
2723: Level: developer
2725: .seealso: DMPlexGetFaceFields()
2726: @*/
2727: PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
2728: {
2732: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);
2733: if (locX_t) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);}
2734: if (locA) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);}
2735: return(0);
2736: }
2738: /*@C
2739: DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces
2741: Input Parameters:
2742: + dm - The DM
2743: . fStart - The first face to include
2744: . fEnd - The first face to exclude
2745: . locX - A local vector with the solution fields
2746: . locX_t - A local vector with solution field time derivatives, or NULL
2747: . faceGeometry - A local vector with face geometry
2748: . cellGeometry - A local vector with cell geometry
2749: - locaGrad - A local vector with field gradients, or NULL
2751: Output Parameters:
2752: + Nface - The number of faces with field values
2753: . uL - The field values at the left side of the face
2754: - uR - The field values at the right side of the face
2756: Level: developer
2758: .seealso: DMPlexGetCellFields()
2759: @*/
2760: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
2761: {
2762: DM dmFace, dmCell, dmGrad = NULL;
2763: PetscSection section;
2764: PetscDS prob;
2765: DMLabel ghostLabel;
2766: const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
2767: PetscBool *isFE;
2768: PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
2769: PetscErrorCode ierr;
2780: DMGetDimension(dm, &dim);
2781: DMGetDS(dm, &prob);
2782: DMGetSection(dm, §ion);
2783: PetscDSGetNumFields(prob, &Nf);
2784: PetscDSGetTotalComponents(prob, &Nc);
2785: PetscMalloc1(Nf, &isFE);
2786: for (f = 0; f < Nf; ++f) {
2787: PetscObject obj;
2788: PetscClassId id;
2790: PetscDSGetDiscretization(prob, f, &obj);
2791: PetscObjectGetClassId(obj, &id);
2792: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
2793: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2794: else {isFE[f] = PETSC_FALSE;}
2795: }
2796: DMGetLabel(dm, "ghost", &ghostLabel);
2797: VecGetArrayRead(locX, &x);
2798: VecGetDM(faceGeometry, &dmFace);
2799: VecGetArrayRead(faceGeometry, &facegeom);
2800: VecGetDM(cellGeometry, &dmCell);
2801: VecGetArrayRead(cellGeometry, &cellgeom);
2802: if (locGrad) {
2803: VecGetDM(locGrad, &dmGrad);
2804: VecGetArrayRead(locGrad, &lgrad);
2805: }
2806: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);
2807: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);
2808: /* Right now just eat the extra work for FE (could make a cell loop) */
2809: for (face = fStart, iface = 0; face < fEnd; ++face) {
2810: const PetscInt *cells;
2811: PetscFVFaceGeom *fg;
2812: PetscFVCellGeom *cgL, *cgR;
2813: PetscScalar *xL, *xR, *gL, *gR;
2814: PetscScalar *uLl = *uL, *uRl = *uR;
2815: PetscInt ghost, nsupp, nchild;
2817: DMLabelGetValue(ghostLabel, face, &ghost);
2818: DMPlexGetSupportSize(dm, face, &nsupp);
2819: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
2820: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
2821: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
2822: DMPlexGetSupport(dm, face, &cells);
2823: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
2824: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
2825: for (f = 0; f < Nf; ++f) {
2826: PetscInt off;
2828: PetscDSGetComponentOffset(prob, f, &off);
2829: if (isFE[f]) {
2830: const PetscInt *cone;
2831: PetscInt comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d;
2833: xL = xR = NULL;
2834: PetscSectionGetFieldComponents(section, f, &comp);
2835: DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
2836: DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
2837: DMPlexGetCone(dm, cells[0], &cone);
2838: DMPlexGetConeSize(dm, cells[0], &coneSizeL);
2839: for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break;
2840: DMPlexGetCone(dm, cells[1], &cone);
2841: DMPlexGetConeSize(dm, cells[1], &coneSizeR);
2842: for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break;
2843: if (faceLocL == coneSizeL && faceLocR == coneSizeR) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d or cell %d", face, cells[0], cells[1]);
2844: /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
2845: /* TODO: this is a hack that might not be right for nonconforming */
2846: if (faceLocL < coneSizeL) {
2847: EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
2848: if (rdof == ldof && faceLocR < coneSizeR) {EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
2849: else {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
2850: }
2851: else {
2852: EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);
2853: PetscSectionGetFieldComponents(section, f, &comp);
2854: for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d];
2855: }
2856: DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
2857: DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
2858: } else {
2859: PetscFV fv;
2860: PetscInt numComp, c;
2862: PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
2863: PetscFVGetNumComponents(fv, &numComp);
2864: DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
2865: DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
2866: if (dmGrad) {
2867: PetscReal dxL[3], dxR[3];
2869: DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
2870: DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
2871: DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
2872: DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
2873: for (c = 0; c < numComp; ++c) {
2874: uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
2875: uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
2876: }
2877: } else {
2878: for (c = 0; c < numComp; ++c) {
2879: uLl[iface*Nc+off+c] = xL[c];
2880: uRl[iface*Nc+off+c] = xR[c];
2881: }
2882: }
2883: }
2884: }
2885: ++iface;
2886: }
2887: *Nface = iface;
2888: VecRestoreArrayRead(locX, &x);
2889: VecRestoreArrayRead(faceGeometry, &facegeom);
2890: VecRestoreArrayRead(cellGeometry, &cellgeom);
2891: if (locGrad) {
2892: VecRestoreArrayRead(locGrad, &lgrad);
2893: }
2894: PetscFree(isFE);
2895: return(0);
2896: }
2898: /*@C
2899: DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces
2901: Input Parameters:
2902: + dm - The DM
2903: . fStart - The first face to include
2904: . fEnd - The first face to exclude
2905: . locX - A local vector with the solution fields
2906: . locX_t - A local vector with solution field time derivatives, or NULL
2907: . faceGeometry - A local vector with face geometry
2908: . cellGeometry - A local vector with cell geometry
2909: - locaGrad - A local vector with field gradients, or NULL
2911: Output Parameters:
2912: + Nface - The number of faces with field values
2913: . uL - The field values at the left side of the face
2914: - uR - The field values at the right side of the face
2916: Level: developer
2918: .seealso: DMPlexGetFaceFields()
2919: @*/
2920: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
2921: {
2925: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);
2926: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);
2927: return(0);
2928: }
2930: /*@C
2931: DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces
2933: Input Parameters:
2934: + dm - The DM
2935: . fStart - The first face to include
2936: . fEnd - The first face to exclude
2937: . faceGeometry - A local vector with face geometry
2938: - cellGeometry - A local vector with cell geometry
2940: Output Parameters:
2941: + Nface - The number of faces with field values
2942: . fgeom - The extract the face centroid and normal
2943: - vol - The cell volume
2945: Level: developer
2947: .seealso: DMPlexGetCellFields()
2948: @*/
2949: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
2950: {
2951: DM dmFace, dmCell;
2952: DMLabel ghostLabel;
2953: const PetscScalar *facegeom, *cellgeom;
2954: PetscInt dim, numFaces = fEnd - fStart, iface, face;
2955: PetscErrorCode ierr;
2963: DMGetDimension(dm, &dim);
2964: DMGetLabel(dm, "ghost", &ghostLabel);
2965: VecGetDM(faceGeometry, &dmFace);
2966: VecGetArrayRead(faceGeometry, &facegeom);
2967: VecGetDM(cellGeometry, &dmCell);
2968: VecGetArrayRead(cellGeometry, &cellgeom);
2969: PetscMalloc1(numFaces, fgeom);
2970: DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);
2971: for (face = fStart, iface = 0; face < fEnd; ++face) {
2972: const PetscInt *cells;
2973: PetscFVFaceGeom *fg;
2974: PetscFVCellGeom *cgL, *cgR;
2975: PetscFVFaceGeom *fgeoml = *fgeom;
2976: PetscReal *voll = *vol;
2977: PetscInt ghost, d, nchild, nsupp;
2979: DMLabelGetValue(ghostLabel, face, &ghost);
2980: DMPlexGetSupportSize(dm, face, &nsupp);
2981: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
2982: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
2983: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
2984: DMPlexGetSupport(dm, face, &cells);
2985: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
2986: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
2987: for (d = 0; d < dim; ++d) {
2988: fgeoml[iface].centroid[d] = fg->centroid[d];
2989: fgeoml[iface].normal[d] = fg->normal[d];
2990: }
2991: voll[iface*2+0] = cgL->volume;
2992: voll[iface*2+1] = cgR->volume;
2993: ++iface;
2994: }
2995: *Nface = iface;
2996: VecRestoreArrayRead(faceGeometry, &facegeom);
2997: VecRestoreArrayRead(cellGeometry, &cellgeom);
2998: return(0);
2999: }
3001: /*@C
3002: DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces
3004: Input Parameters:
3005: + dm - The DM
3006: . fStart - The first face to include
3007: . fEnd - The first face to exclude
3008: . faceGeometry - A local vector with face geometry
3009: - cellGeometry - A local vector with cell geometry
3011: Output Parameters:
3012: + Nface - The number of faces with field values
3013: . fgeom - The extract the face centroid and normal
3014: - vol - The cell volume
3016: Level: developer
3018: .seealso: DMPlexGetFaceFields()
3019: @*/
3020: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3021: {
3025: PetscFree(*fgeom);
3026: DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);
3027: return(0);
3028: }
3030: PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3031: {
3032: char composeStr[33] = {0};
3033: PetscObjectId id;
3034: PetscContainer container;
3035: PetscErrorCode ierr;
3038: PetscObjectGetId((PetscObject)quad,&id);
3039: PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);
3040: PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
3041: if (container) {
3042: PetscContainerGetPointer(container, (void **) geom);
3043: } else {
3044: DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
3045: PetscContainerCreate(PETSC_COMM_SELF,&container);
3046: PetscContainerSetPointer(container, (void *) *geom);
3047: PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
3048: PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
3049: PetscContainerDestroy(&container);
3050: }
3051: return(0);
3052: }
3054: PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3055: {
3057: *geom = NULL;
3058: return(0);
3059: }
3061: PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user)
3062: {
3063: DM_Plex *mesh = (DM_Plex *) dm->data;
3064: const char *name = "Residual";
3065: DM dmAux = NULL;
3066: DMLabel ghostLabel = NULL;
3067: PetscDS prob = NULL;
3068: PetscDS probAux = NULL;
3069: PetscBool useFEM = PETSC_FALSE;
3070: PetscBool isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
3071: DMField coordField = NULL;
3072: Vec locA;
3073: PetscScalar *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL;
3074: IS chunkIS;
3075: const PetscInt *cells;
3076: PetscInt cStart, cEnd, numCells;
3077: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd;
3078: PetscInt maxDegree = PETSC_MAX_INT;
3079: PetscQuadrature affineQuad = NULL, *quads = NULL;
3080: PetscFEGeom *affineGeom = NULL, **geoms = NULL;
3081: PetscErrorCode ierr;
3084: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
3085: /* FEM+FVM */
3086: /* 1: Get sizes from dm and dmAux */
3087: DMGetLabel(dm, "ghost", &ghostLabel);
3088: DMGetDS(dm, &prob);
3089: PetscDSGetNumFields(prob, &Nf);
3090: PetscDSGetTotalDimension(prob, &totDim);
3091: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
3092: if (locA) {
3093: VecGetDM(locA, &dmAux);
3094: DMGetDS(dmAux, &probAux);
3095: PetscDSGetTotalDimension(probAux, &totDimAux);
3096: }
3097: /* 2: Get geometric data */
3098: for (f = 0; f < Nf; ++f) {
3099: PetscObject obj;
3100: PetscClassId id;
3101: PetscBool fimp;
3103: PetscDSGetImplicit(prob, f, &fimp);
3104: if (isImplicit != fimp) continue;
3105: PetscDSGetDiscretization(prob, f, &obj);
3106: PetscObjectGetClassId(obj, &id);
3107: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
3108: if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented");
3109: }
3110: if (useFEM) {
3111: DMGetCoordinateField(dm, &coordField);
3112: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
3113: if (maxDegree <= 1) {
3114: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
3115: if (affineQuad) {
3116: DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3117: }
3118: } else {
3119: PetscCalloc2(Nf,&quads,Nf,&geoms);
3120: for (f = 0; f < Nf; ++f) {
3121: PetscObject obj;
3122: PetscClassId id;
3123: PetscBool fimp;
3125: PetscDSGetImplicit(prob, f, &fimp);
3126: if (isImplicit != fimp) continue;
3127: PetscDSGetDiscretization(prob, f, &obj);
3128: PetscObjectGetClassId(obj, &id);
3129: if (id == PETSCFE_CLASSID) {
3130: PetscFE fe = (PetscFE) obj;
3132: PetscFEGetQuadrature(fe, &quads[f]);
3133: PetscObjectReference((PetscObject)quads[f]);
3134: DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3135: }
3136: }
3137: }
3138: }
3139: /* Loop over chunks */
3140: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3141: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
3142: if (useFEM) {ISCreate(PETSC_COMM_SELF, &chunkIS);}
3143: numCells = cEnd - cStart;
3144: numChunks = 1;
3145: cellChunkSize = numCells/numChunks;
3146: numChunks = PetscMin(1,numCells);
3147: for (chunk = 0; chunk < numChunks; ++chunk) {
3148: PetscScalar *elemVec, *fluxL = NULL, *fluxR = NULL;
3149: PetscReal *vol = NULL;
3150: PetscFVFaceGeom *fgeom = NULL;
3151: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
3152: PetscInt numFaces = 0;
3154: /* Extract field coefficients */
3155: if (useFEM) {
3156: ISGetPointSubrange(chunkIS, cS, cE, cells);
3157: DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3158: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3159: PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
3160: }
3161: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
3162: /* Loop over fields */
3163: for (f = 0; f < Nf; ++f) {
3164: PetscObject obj;
3165: PetscClassId id;
3166: PetscBool fimp;
3167: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
3169: PetscDSGetImplicit(prob, f, &fimp);
3170: if (isImplicit != fimp) continue;
3171: PetscDSGetDiscretization(prob, f, &obj);
3172: PetscObjectGetClassId(obj, &id);
3173: if (id == PETSCFE_CLASSID) {
3174: PetscFE fe = (PetscFE) obj;
3175: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
3176: PetscFEGeom *chunkGeom = NULL;
3177: PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
3178: PetscInt Nq, Nb;
3180: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3181: PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
3182: PetscFEGetDimension(fe, &Nb);
3183: blockSize = Nb;
3184: batchSize = numBlocks * blockSize;
3185: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3186: numChunks = numCells / (numBatches*batchSize);
3187: Ne = numChunks*numBatches*batchSize;
3188: Nr = numCells % (numBatches*batchSize);
3189: offset = numCells - Nr;
3190: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
3191: /* For FV, I think we use a P0 basis and the cell coefficients (for subdivided cells, we can tweak the basis tabulation to be the indicator function) */
3192: PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
3193: PetscFEIntegrateResidual(fe, prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
3194: PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
3195: PetscFEIntegrateResidual(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
3196: PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
3197: } else if (id == PETSCFV_CLASSID) {
3198: PetscFV fv = (PetscFV) obj;
3200: Ne = numFaces;
3201: /* Riemann solve over faces (need fields at face centroids) */
3202: /* We need to evaluate FE fields at those coordinates */
3203: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
3204: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
3205: }
3206: /* Loop over domain */
3207: if (useFEM) {
3208: /* Add elemVec to locX */
3209: for (c = cS; c < cE; ++c) {
3210: const PetscInt cell = cells ? cells[c] : c;
3211: const PetscInt cind = c - cStart;
3213: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);}
3214: if (ghostLabel) {
3215: PetscInt ghostVal;
3217: DMLabelGetValue(ghostLabel,cell,&ghostVal);
3218: if (ghostVal > 0) continue;
3219: }
3220: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
3221: }
3222: }
3223: /* Handle time derivative */
3224: if (locX_t) {
3225: PetscScalar *x_t, *fa;
3227: VecGetArray(locF, &fa);
3228: VecGetArray(locX_t, &x_t);
3229: for (f = 0; f < Nf; ++f) {
3230: PetscFV fv;
3231: PetscObject obj;
3232: PetscClassId id;
3233: PetscInt pdim, d;
3235: PetscDSGetDiscretization(prob, f, &obj);
3236: PetscObjectGetClassId(obj, &id);
3237: if (id != PETSCFV_CLASSID) continue;
3238: fv = (PetscFV) obj;
3239: PetscFVGetNumComponents(fv, &pdim);
3240: for (c = cS; c < cE; ++c) {
3241: const PetscInt cell = cells ? cells[c] : c;
3242: PetscScalar *u_t, *r;
3244: if (ghostLabel) {
3245: PetscInt ghostVal;
3247: DMLabelGetValue(ghostLabel, cell, &ghostVal);
3248: if (ghostVal > 0) continue;
3249: }
3250: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
3251: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
3252: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
3253: }
3254: }
3255: VecRestoreArray(locX_t, &x_t);
3256: VecRestoreArray(locF, &fa);
3257: }
3258: if (useFEM) {
3259: DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3260: DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3261: }
3262: }
3263: if (useFEM) {ISDestroy(&chunkIS);}
3264: ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3265: /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */
3266: if (useFEM) {
3267: if (maxDegree <= 1) {
3268: DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3269: PetscQuadratureDestroy(&affineQuad);
3270: } else {
3271: for (f = 0; f < Nf; ++f) {
3272: DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3273: PetscQuadratureDestroy(&quads[f]);
3274: }
3275: PetscFree2(quads,geoms);
3276: }
3277: }
3278: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
3279: return(0);
3280: }
3282: /*
3283: We always assemble JacP, and if the matrix is different from Jac and two different sets of point functions are provided, we also assemble Jac
3285: X - The local solution vector
3286: X_t - The local solution time derviative vector, or NULL
3287: */
3288: PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS,
3289: PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx)
3290: {
3291: DM_Plex *mesh = (DM_Plex *) dm->data;
3292: const char *name = "Jacobian", *nameP = "JacobianPre";
3293: DM dmAux = NULL;
3294: PetscDS prob, probAux = NULL;
3295: PetscSection sectionAux = NULL;
3296: Vec A;
3297: DMField coordField;
3298: PetscFEGeom *cgeomFEM;
3299: PetscQuadrature qGeom = NULL;
3300: Mat J = Jac, JP = JacP;
3301: PetscScalar *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL;
3302: PetscBool hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE;
3303: const PetscInt *cells;
3304: PetscInt Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0;
3305: PetscErrorCode ierr;
3308: CHKMEMQ;
3309: ISGetLocalSize(cellIS, &numCells);
3310: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3311: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
3312: DMGetDS(dm, &prob);
3313: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3314: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3315: if (dmAux) {
3316: DMGetDefaultSection(dmAux, §ionAux);
3317: DMGetDS(dmAux, &probAux);
3318: }
3319: /* Get flags */
3320: PetscDSGetNumFields(prob, &Nf);
3321: DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3322: for (fieldI = 0; fieldI < Nf; ++fieldI) {
3323: PetscObject disc;
3324: PetscClassId id;
3325: PetscDSGetDiscretization(prob, fieldI, &disc);
3326: PetscObjectGetClassId(disc, &id);
3327: if (id == PETSCFE_CLASSID) {isFE[fieldI] = PETSC_TRUE;}
3328: else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;}
3329: }
3330: PetscDSHasJacobian(prob, &hasJac);
3331: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
3332: PetscDSHasDynamicJacobian(prob, &hasDyn);
3333: assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE;
3334: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
3335: PetscObjectTypeCompare((PetscObject) Jac, MATIS, &isMatIS);
3336: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
3337: /* Setup input data and temp arrays (should be DMGetWorkArray) */
3338: if (isMatISP || isMatISP) {DMPlexGetSubdomainSection(dm, &globalSection);}
3339: if (isMatIS) {MatISGetLocalMat(Jac, &J);}
3340: if (isMatISP) {MatISGetLocalMat(JacP, &JP);}
3341: if (hasFV) {MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);} /* No allocated space for FV stuff, so ignore the zero entries */
3342: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3343: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3344: PetscDSGetTotalDimension(prob, &totDim);
3345: if (probAux) {PetscDSGetTotalDimension(probAux, &totDimAux);}
3346: CHKMEMQ;
3347: /* Compute batch sizes */
3348: if (isFE[0]) {
3349: PetscFE fe;
3350: PetscQuadrature q;
3351: PetscInt numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb;
3353: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3354: PetscFEGetQuadrature(fe, &q);
3355: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
3356: PetscFEGetDimension(fe, &Nb);
3357: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3358: blockSize = Nb*numQuadPoints;
3359: batchSize = numBlocks * blockSize;
3360: chunkSize = numBatches * batchSize;
3361: numChunks = numCells / chunkSize + numCells % chunkSize;
3362: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3363: } else {
3364: chunkSize = numCells;
3365: numChunks = 1;
3366: }
3367: /* Get work space */
3368: wsz = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize;
3369: DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);
3370: PetscMemzero(work, wsz * sizeof(PetscScalar));
3371: off = 0;
3372: u = X ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
3373: u_t = X_t ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
3374: a = dmAux ? (sz = chunkSize*totDimAux, off += sz, work+off-sz) : NULL;
3375: elemMat = hasJac ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3376: elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3377: elemMatD = hasDyn ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3378: if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz);
3379: /* Setup geometry */
3380: DMGetCoordinateField(dm, &coordField);
3381: DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
3382: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);}
3383: if (!qGeom) {
3384: PetscFE fe;
3386: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3387: PetscFEGetQuadrature(fe, &qGeom);
3388: PetscObjectReference((PetscObject) qGeom);
3389: }
3390: DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3391: /* Compute volume integrals */
3392: if (assembleJac) {MatZeroEntries(J);}
3393: MatZeroEntries(JP);
3394: for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) {
3395: const PetscInt Ncell = PetscMin(chunkSize, numCells - offCell);
3396: PetscInt c;
3398: /* Extract values */
3399: for (c = 0; c < Ncell; ++c) {
3400: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3401: PetscScalar *x = NULL, *x_t = NULL;
3402: PetscInt i;
3404: if (X) {
3405: DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
3406: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
3407: DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
3408: }
3409: if (X_t) {
3410: DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
3411: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
3412: DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
3413: }
3414: if (dmAux) {
3415: DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);
3416: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
3417: DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);
3418: }
3419: }
3420: CHKMEMQ;
3421: for (fieldI = 0; fieldI < Nf; ++fieldI) {
3422: PetscFE fe;
3423: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3424: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
3425: if (hasJac) {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);}
3426: if (hasPrec) {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);}
3427: if (hasDyn) {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);}
3428: }
3429: /* For finite volume, add the identity */
3430: if (!isFE[fieldI]) {
3431: PetscFV fv;
3432: PetscInt eOffset = 0, Nc, fc, foff;
3434: PetscDSGetFieldOffset(prob, fieldI, &foff);
3435: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
3436: PetscFVGetNumComponents(fv, &Nc);
3437: for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) {
3438: for (fc = 0; fc < Nc; ++fc) {
3439: const PetscInt i = foff + fc;
3440: if (hasJac) {elemMat [eOffset+i*totDim+i] = 1.0;}
3441: if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;}
3442: }
3443: }
3444: }
3445: }
3446: CHKMEMQ;
3447: /* Add contribution from X_t */
3448: if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
3449: /* Insert values into matrix */
3450: for (c = 0; c < Ncell; ++c) {
3451: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3452: if (mesh->printFEM > 1) {
3453: if (hasJac) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
3454: if (hasPrec) {DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
3455: }
3456: if (assembleJac) {DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);}
3457: DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
3458: }
3459: CHKMEMQ;
3460: }
3461: /* Cleanup */
3462: DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3463: PetscQuadratureDestroy(&qGeom);
3464: if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
3465: DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3466: DMRestoreWorkArray(dm, ((1 + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize, MPIU_SCALAR, &work);
3467: /* Compute boundary integrals */
3468: /* DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx); */
3469: /* Assemble matrix */
3470: if (assembleJac) {MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);}
3471: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
3472: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
3473: CHKMEMQ;
3474: return(0);
3475: }