Actual source code: plexfem.c
petsc-3.10.5 2019-03-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 PetscContainerUserDestroy_PetscFEGeom (void *ctx)
9: {
10: PetscFEGeom *geom = (PetscFEGeom *) ctx;
14: PetscFEGeomDestroy(&geom);
15: return(0);
16: }
18: static PetscErrorCode DMPlexGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
19: {
20: char composeStr[33] = {0};
21: PetscObjectId id;
22: PetscContainer container;
23: PetscErrorCode ierr;
26: PetscObjectGetId((PetscObject)quad,&id);
27: PetscSNPrintf(composeStr, 32, "DMPlexGetFEGeom_%x\n", id);
28: PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
29: if (container) {
30: PetscContainerGetPointer(container, (void **) geom);
31: } else {
32: DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
33: PetscContainerCreate(PETSC_COMM_SELF,&container);
34: PetscContainerSetPointer(container, (void *) *geom);
35: PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
36: PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
37: PetscContainerDestroy(&container);
38: }
39: return(0);
40: }
42: static PetscErrorCode DMPlexRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
43: {
45: *geom = NULL;
46: return(0);
47: }
49: /*@
50: DMPlexGetScale - Get the scale for the specified fundamental unit
52: Not collective
54: Input Arguments:
55: + dm - the DM
56: - unit - The SI unit
58: Output Argument:
59: . scale - The value used to scale all quantities with this unit
61: Level: advanced
63: .seealso: DMPlexSetScale(), PetscUnit
64: @*/
65: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
66: {
67: DM_Plex *mesh = (DM_Plex*) dm->data;
72: *scale = mesh->scale[unit];
73: return(0);
74: }
76: /*@
77: DMPlexSetScale - Set the scale for the specified fundamental unit
79: Not collective
81: Input Arguments:
82: + dm - the DM
83: . unit - The SI unit
84: - scale - The value used to scale all quantities with this unit
86: Level: advanced
88: .seealso: DMPlexGetScale(), PetscUnit
89: @*/
90: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
91: {
92: DM_Plex *mesh = (DM_Plex*) dm->data;
96: mesh->scale[unit] = scale;
97: return(0);
98: }
100: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
101: {
102: 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}}};
103: PetscInt *ctxInt = (PetscInt *) ctx;
104: PetscInt dim2 = ctxInt[0];
105: PetscInt d = ctxInt[1];
106: PetscInt i, j, k = dim > 2 ? d - dim : d;
109: if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
110: for (i = 0; i < dim; i++) mode[i] = 0.;
111: if (d < dim) {
112: mode[d] = 1.; /* Translation along axis d */
113: } else {
114: for (i = 0; i < dim; i++) {
115: for (j = 0; j < dim; j++) {
116: mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */
117: }
118: }
119: }
120: return(0);
121: }
123: /*@C
124: DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation
126: Collective on DM
128: Input Arguments:
129: . dm - the DM
131: Output Argument:
132: . sp - the null space
134: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
136: Level: advanced
138: .seealso: MatNullSpaceCreate(), PCGAMG
139: @*/
140: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
141: {
142: MPI_Comm comm;
143: Vec mode[6];
144: PetscSection section, globalSection;
145: PetscInt dim, dimEmbed, n, m, d, i, j;
149: PetscObjectGetComm((PetscObject)dm,&comm);
150: DMGetDimension(dm, &dim);
151: DMGetCoordinateDim(dm, &dimEmbed);
152: if (dim == 1) {
153: MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
154: return(0);
155: }
156: DMGetSection(dm, §ion);
157: DMGetGlobalSection(dm, &globalSection);
158: PetscSectionGetConstrainedStorageSize(globalSection, &n);
159: m = (dim*(dim+1))/2;
160: VecCreate(comm, &mode[0]);
161: VecSetSizes(mode[0], n, PETSC_DETERMINE);
162: VecSetUp(mode[0]);
163: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
164: for (d = 0; d < m; d++) {
165: PetscInt ctx[2];
166: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
167: void *voidctx = (void *) (&ctx[0]);
169: ctx[0] = dimEmbed;
170: ctx[1] = d;
171: DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
172: }
173: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
174: /* Orthonormalize system */
175: for (i = dim; i < m; ++i) {
176: PetscScalar dots[6];
178: VecMDot(mode[i], i, mode, dots);
179: for (j = 0; j < i; ++j) dots[j] *= -1.0;
180: VecMAXPY(mode[i], i, dots, mode);
181: VecNormalize(mode[i], NULL);
182: }
183: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
184: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
185: return(0);
186: }
188: /*@C
189: DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation
191: Collective on DM
193: Input Arguments:
194: + dm - the DM
195: . nb - The number of bodies
196: . label - The DMLabel marking each domain
197: . nids - The number of ids per body
198: - ids - An array of the label ids in sequence for each domain
200: Output Argument:
201: . sp - the null space
203: Note: This is necessary to provide a suitable coarse space for algebraic multigrid
205: Level: advanced
207: .seealso: MatNullSpaceCreate()
208: @*/
209: PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp)
210: {
211: MPI_Comm comm;
212: PetscSection section, globalSection;
213: Vec *mode;
214: PetscScalar *dots;
215: PetscInt dim, dimEmbed, n, m, b, d, i, j, off;
219: PetscObjectGetComm((PetscObject)dm,&comm);
220: DMGetDimension(dm, &dim);
221: DMGetCoordinateDim(dm, &dimEmbed);
222: DMGetSection(dm, §ion);
223: DMGetGlobalSection(dm, &globalSection);
224: PetscSectionGetConstrainedStorageSize(globalSection, &n);
225: m = nb * (dim*(dim+1))/2;
226: PetscMalloc2(m, &mode, m, &dots);
227: VecCreate(comm, &mode[0]);
228: VecSetSizes(mode[0], n, PETSC_DETERMINE);
229: VecSetUp(mode[0]);
230: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
231: for (b = 0, off = 0; b < nb; ++b) {
232: for (d = 0; d < m/nb; ++d) {
233: PetscInt ctx[2];
234: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
235: void *voidctx = (void *) (&ctx[0]);
237: ctx[0] = dimEmbed;
238: ctx[1] = d;
239: DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
240: off += nids[b];
241: }
242: }
243: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
244: /* Orthonormalize system */
245: for (i = 0; i < m; ++i) {
246: VecMDot(mode[i], i, mode, dots);
247: for (j = 0; j < i; ++j) dots[j] *= -1.0;
248: VecMAXPY(mode[i], i, dots, mode);
249: VecNormalize(mode[i], NULL);
250: }
251: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
252: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
253: PetscFree2(mode, dots);
254: return(0);
255: }
257: /*@
258: DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
259: are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
260: evaluating the dual space basis of that point. A basis function is associated with the point in its
261: transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
262: projection height, which is set with this function. By default, the maximum projection height is zero, which means
263: that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior
264: basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
266: Input Parameters:
267: + dm - the DMPlex object
268: - height - the maximum projection height >= 0
270: Level: advanced
272: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
273: @*/
274: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
275: {
276: DM_Plex *plex = (DM_Plex *) dm->data;
280: plex->maxProjectionHeight = height;
281: return(0);
282: }
284: /*@
285: DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
286: DMPlexProjectXXXLocal() functions.
288: Input Parameters:
289: . dm - the DMPlex object
291: Output Parameters:
292: . height - the maximum projection height
294: Level: intermediate
296: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
297: @*/
298: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
299: {
300: DM_Plex *plex = (DM_Plex *) dm->data;
304: *height = plex->maxProjectionHeight;
305: return(0);
306: }
308: /*@C
309: DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector
311: Input Parameters:
312: + dm - The DM, with a PetscDS that matches the problem being constrained
313: . time - The time
314: . field - The field to constrain
315: . Nc - The number of constrained field components, or 0 for all components
316: . comps - An array of constrained component numbers, or NULL for all components
317: . label - The DMLabel defining constrained points
318: . numids - The number of DMLabel ids for constrained points
319: . ids - An array of ids for constrained points
320: . func - A pointwise function giving boundary values
321: - ctx - An optional user context for bcFunc
323: Output Parameter:
324: . locX - A local vector to receives the boundary values
326: Level: developer
328: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
329: @*/
330: 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)
331: {
332: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
333: void **ctxs;
334: PetscInt numFields;
335: PetscErrorCode ierr;
338: DMGetNumFields(dm, &numFields);
339: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
340: funcs[field] = func;
341: ctxs[field] = ctx;
342: DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
343: PetscFree2(funcs,ctxs);
344: return(0);
345: }
347: /*@C
348: DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector
350: Input Parameters:
351: + dm - The DM, with a PetscDS that matches the problem being constrained
352: . time - The time
353: . locU - A local vector with the input solution values
354: . field - The field to constrain
355: . Nc - The number of constrained field components, or 0 for all components
356: . comps - An array of constrained component numbers, or NULL for all components
357: . label - The DMLabel defining constrained points
358: . numids - The number of DMLabel ids for constrained points
359: . ids - An array of ids for constrained points
360: . func - A pointwise function giving boundary values
361: - ctx - An optional user context for bcFunc
363: Output Parameter:
364: . locX - A local vector to receives the boundary values
366: Level: developer
368: .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary()
369: @*/
370: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
371: void (*func)(PetscInt, PetscInt, PetscInt,
372: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
373: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
374: PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
375: PetscScalar[]),
376: void *ctx, Vec locX)
377: {
378: void (**funcs)(PetscInt, PetscInt, PetscInt,
379: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
380: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
381: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
382: void **ctxs;
383: PetscInt numFields;
384: PetscErrorCode ierr;
387: DMGetNumFields(dm, &numFields);
388: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
389: funcs[field] = func;
390: ctxs[field] = ctx;
391: DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
392: PetscFree2(funcs,ctxs);
393: return(0);
394: }
396: /*@C
397: DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector
399: Input Parameters:
400: + dm - The DM, with a PetscDS that matches the problem being constrained
401: . time - The time
402: . faceGeometry - A vector with the FVM face geometry information
403: . cellGeometry - A vector with the FVM cell geometry information
404: . Grad - A vector with the FVM cell gradient information
405: . field - The field to constrain
406: . Nc - The number of constrained field components, or 0 for all components
407: . comps - An array of constrained component numbers, or NULL for all components
408: . label - The DMLabel defining constrained points
409: . numids - The number of DMLabel ids for constrained points
410: . ids - An array of ids for constrained points
411: . func - A pointwise function giving boundary values
412: - ctx - An optional user context for bcFunc
414: Output Parameter:
415: . locX - A local vector to receives the boundary values
417: Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary()
419: Level: developer
421: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
422: @*/
423: 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[],
424: PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
425: {
426: PetscDS prob;
427: PetscSF sf;
428: DM dmFace, dmCell, dmGrad;
429: const PetscScalar *facegeom, *cellgeom = NULL, *grad;
430: const PetscInt *leaves;
431: PetscScalar *x, *fx;
432: PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i;
433: PetscErrorCode ierr, ierru = 0;
436: DMGetPointSF(dm, &sf);
437: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
438: nleaves = PetscMax(0, nleaves);
439: DMGetDimension(dm, &dim);
440: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
441: DMGetDS(dm, &prob);
442: VecGetDM(faceGeometry, &dmFace);
443: VecGetArrayRead(faceGeometry, &facegeom);
444: if (cellGeometry) {
445: VecGetDM(cellGeometry, &dmCell);
446: VecGetArrayRead(cellGeometry, &cellgeom);
447: }
448: if (Grad) {
449: PetscFV fv;
451: PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
452: VecGetDM(Grad, &dmGrad);
453: VecGetArrayRead(Grad, &grad);
454: PetscFVGetNumComponents(fv, &pdim);
455: DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
456: }
457: VecGetArray(locX, &x);
458: for (i = 0; i < numids; ++i) {
459: IS faceIS;
460: const PetscInt *faces;
461: PetscInt numFaces, f;
463: DMLabelGetStratumIS(label, ids[i], &faceIS);
464: if (!faceIS) continue; /* No points with that id on this process */
465: ISGetLocalSize(faceIS, &numFaces);
466: ISGetIndices(faceIS, &faces);
467: for (f = 0; f < numFaces; ++f) {
468: const PetscInt face = faces[f], *cells;
469: PetscFVFaceGeom *fg;
471: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
472: PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
473: if (loc >= 0) continue;
474: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
475: DMPlexGetSupport(dm, face, &cells);
476: if (Grad) {
477: PetscFVCellGeom *cg;
478: PetscScalar *cx, *cgrad;
479: PetscScalar *xG;
480: PetscReal dx[3];
481: PetscInt d;
483: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
484: DMPlexPointLocalRead(dm, cells[0], x, &cx);
485: DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
486: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
487: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
488: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
489: ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
490: if (ierru) {
491: ISRestoreIndices(faceIS, &faces);
492: ISDestroy(&faceIS);
493: goto cleanup;
494: }
495: } else {
496: PetscScalar *xI;
497: PetscScalar *xG;
499: DMPlexPointLocalRead(dm, cells[0], x, &xI);
500: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
501: ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
502: if (ierru) {
503: ISRestoreIndices(faceIS, &faces);
504: ISDestroy(&faceIS);
505: goto cleanup;
506: }
507: }
508: }
509: ISRestoreIndices(faceIS, &faces);
510: ISDestroy(&faceIS);
511: }
512: cleanup:
513: VecRestoreArray(locX, &x);
514: if (Grad) {
515: DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
516: VecRestoreArrayRead(Grad, &grad);
517: }
518: if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
519: VecRestoreArrayRead(faceGeometry, &facegeom);
520: CHKERRQ(ierru);
521: return(0);
522: }
524: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
525: {
526: PetscInt numBd, b;
530: PetscDSGetNumBoundary(dm->prob, &numBd);
531: for (b = 0; b < numBd; ++b) {
532: DMBoundaryConditionType type;
533: const char *labelname;
534: DMLabel label;
535: PetscInt field, Nc;
536: const PetscInt *comps;
537: PetscObject obj;
538: PetscClassId id;
539: void (*func)(void);
540: PetscInt numids;
541: const PetscInt *ids;
542: void *ctx;
544: DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
545: if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
546: DMGetLabel(dm, labelname, &label);
547: DMGetField(dm, field, &obj);
548: PetscObjectGetClassId(obj, &id);
549: if (id == PETSCFE_CLASSID) {
550: switch (type) {
551: /* for FEM, there is no insertion to be done for non-essential boundary conditions */
552: case DM_BC_ESSENTIAL:
553: DMPlexLabelAddCells(dm,label);
554: DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
555: DMPlexLabelClearCells(dm,label);
556: break;
557: case DM_BC_ESSENTIAL_FIELD:
558: DMPlexLabelAddCells(dm,label);
559: DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
560: (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
561: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
562: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
563: DMPlexLabelClearCells(dm,label);
564: break;
565: default: break;
566: }
567: } else if (id == PETSCFV_CLASSID) {
568: if (!faceGeomFVM) continue;
569: DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
570: (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
571: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
572: }
573: return(0);
574: }
576: /*@
577: DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector
579: Input Parameters:
580: + dm - The DM
581: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
582: . time - The time
583: . faceGeomFVM - Face geometry data for FV discretizations
584: . cellGeomFVM - Cell geometry data for FV discretizations
585: - gradFVM - Gradient reconstruction data for FV discretizations
587: Output Parameters:
588: . locX - Solution updated with boundary values
590: Level: developer
592: .seealso: DMProjectFunctionLabelLocal()
593: @*/
594: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
595: {
604: PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
605: return(0);
606: }
608: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
609: {
610: Vec localX;
611: PetscErrorCode ierr;
614: DMGetLocalVector(dm, &localX);
615: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
616: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
617: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
618: DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
619: DMRestoreLocalVector(dm, &localX);
620: return(0);
621: }
623: /*@C
624: DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
626: Input Parameters:
627: + dm - The DM
628: . time - The time
629: . funcs - The functions to evaluate for each field component
630: . ctxs - Optional array of contexts to pass to each function, or NULL.
631: - localX - The coefficient vector u_h, a local vector
633: Output Parameter:
634: . diff - The diff ||u - u_h||_2
636: Level: developer
638: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
639: @*/
640: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
641: {
642: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
643: PetscSection section;
644: PetscQuadrature quad;
645: PetscScalar *funcVal, *interpolant;
646: PetscReal *coords, *detJ, *J;
647: PetscReal localDiff = 0.0;
648: const PetscReal *quadWeights;
649: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset;
650: PetscErrorCode ierr;
653: DMGetDimension(dm, &dim);
654: DMGetCoordinateDim(dm, &coordDim);
655: DMGetSection(dm, §ion);
656: PetscSectionGetNumFields(section, &numFields);
657: for (field = 0; field < numFields; ++field) {
658: PetscObject obj;
659: PetscClassId id;
660: PetscInt Nc;
662: DMGetField(dm, field, &obj);
663: PetscObjectGetClassId(obj, &id);
664: if (id == PETSCFE_CLASSID) {
665: PetscFE fe = (PetscFE) obj;
667: PetscFEGetQuadrature(fe, &quad);
668: PetscFEGetNumComponents(fe, &Nc);
669: } else if (id == PETSCFV_CLASSID) {
670: PetscFV fv = (PetscFV) obj;
672: PetscFVGetQuadrature(fv, &quad);
673: PetscFVGetNumComponents(fv, &Nc);
674: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
675: numComponents += Nc;
676: }
677: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
678: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
679: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
680: DMPlexGetVTKCellHeight(dm, &cellHeight);
681: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
682: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
683: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
684: for (c = cStart; c < cEnd; ++c) {
685: PetscScalar *x = NULL;
686: PetscReal elemDiff = 0.0;
687: PetscInt qc = 0;
689: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
690: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
692: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
693: PetscObject obj;
694: PetscClassId id;
695: void * const ctx = ctxs ? ctxs[field] : NULL;
696: PetscInt Nb, Nc, q, fc;
698: DMGetField(dm, field, &obj);
699: PetscObjectGetClassId(obj, &id);
700: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
701: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
702: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
703: if (debug) {
704: char title[1024];
705: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
706: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
707: }
708: for (q = 0; q < Nq; ++q) {
709: 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);
710: (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
711: if (ierr) {
712: PetscErrorCode ierr2;
713: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
714: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
715: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
716:
717: }
718: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
719: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
720: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
721: for (fc = 0; fc < Nc; ++fc) {
722: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
723: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
724: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
725: }
726: }
727: fieldOffset += Nb;
728: qc += Nc;
729: }
730: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
731: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
732: localDiff += elemDiff;
733: }
734: PetscFree5(funcVal,interpolant,coords,detJ,J);
735: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
736: *diff = PetscSqrtReal(*diff);
737: return(0);
738: }
740: 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)
741: {
742: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
743: PetscSection section;
744: PetscQuadrature quad;
745: Vec localX;
746: PetscScalar *funcVal, *interpolant;
747: const PetscReal *quadPoints, *quadWeights;
748: PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ;
749: PetscReal localDiff = 0.0;
750: PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
751: PetscErrorCode ierr;
754: DMGetDimension(dm, &dim);
755: DMGetCoordinateDim(dm, &coordDim);
756: DMGetSection(dm, §ion);
757: PetscSectionGetNumFields(section, &numFields);
758: DMGetLocalVector(dm, &localX);
759: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
760: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
761: for (field = 0; field < numFields; ++field) {
762: PetscFE fe;
763: PetscInt Nc;
765: DMGetField(dm, field, (PetscObject *) &fe);
766: PetscFEGetQuadrature(fe, &quad);
767: PetscFEGetNumComponents(fe, &Nc);
768: numComponents += Nc;
769: }
770: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
771: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
772: /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
773: PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
774: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
775: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
776: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
777: for (c = cStart; c < cEnd; ++c) {
778: PetscScalar *x = NULL;
779: PetscReal elemDiff = 0.0;
780: PetscInt qc = 0;
782: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);
783: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
785: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
786: PetscFE fe;
787: void * const ctx = ctxs ? ctxs[field] : NULL;
788: PetscReal *basisDer;
789: PetscInt Nb, Nc, q, fc;
791: DMGetField(dm, field, (PetscObject *) &fe);
792: PetscFEGetDimension(fe, &Nb);
793: PetscFEGetNumComponents(fe, &Nc);
794: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
795: if (debug) {
796: char title[1024];
797: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
798: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
799: }
800: for (q = 0; q < Nq; ++q) {
801: 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);
802: (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
803: if (ierr) {
804: PetscErrorCode ierr2;
805: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
806: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
807: ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
808:
809: }
810: PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
811: for (fc = 0; fc < Nc; ++fc) {
812: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
813: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
814: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
815: }
816: }
817: fieldOffset += Nb;
818: qc += Nc;
819: }
820: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
821: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
822: localDiff += elemDiff;
823: }
824: PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
825: DMRestoreLocalVector(dm, &localX);
826: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
827: *diff = PetscSqrtReal(*diff);
828: return(0);
829: }
831: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
832: {
833: const PetscInt debug = ((DM_Plex*)dm->data)->printL2;
834: PetscSection section;
835: PetscQuadrature quad;
836: Vec localX;
837: PetscScalar *funcVal, *interpolant;
838: PetscReal *coords, *detJ, *J;
839: PetscReal *localDiff;
840: const PetscReal *quadPoints, *quadWeights;
841: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
842: PetscErrorCode ierr;
845: DMGetDimension(dm, &dim);
846: DMGetCoordinateDim(dm, &coordDim);
847: DMGetSection(dm, §ion);
848: PetscSectionGetNumFields(section, &numFields);
849: DMGetLocalVector(dm, &localX);
850: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
851: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
852: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
853: for (field = 0; field < numFields; ++field) {
854: PetscObject obj;
855: PetscClassId id;
856: PetscInt Nc;
858: DMGetField(dm, field, &obj);
859: PetscObjectGetClassId(obj, &id);
860: if (id == PETSCFE_CLASSID) {
861: PetscFE fe = (PetscFE) obj;
863: PetscFEGetQuadrature(fe, &quad);
864: PetscFEGetNumComponents(fe, &Nc);
865: } else if (id == PETSCFV_CLASSID) {
866: PetscFV fv = (PetscFV) obj;
868: PetscFVGetQuadrature(fv, &quad);
869: PetscFVGetNumComponents(fv, &Nc);
870: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
871: numComponents += Nc;
872: }
873: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
874: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
875: PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
876: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
877: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
878: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
879: for (c = cStart; c < cEnd; ++c) {
880: PetscScalar *x = NULL;
881: PetscInt qc = 0;
883: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
884: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
886: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
887: PetscObject obj;
888: PetscClassId id;
889: void * const ctx = ctxs ? ctxs[field] : NULL;
890: PetscInt Nb, Nc, q, fc;
892: PetscReal elemDiff = 0.0;
894: DMGetField(dm, field, &obj);
895: PetscObjectGetClassId(obj, &id);
896: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
897: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
898: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
899: if (debug) {
900: char title[1024];
901: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
902: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
903: }
904: for (q = 0; q < Nq; ++q) {
905: 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);
906: (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
907: if (ierr) {
908: PetscErrorCode ierr2;
909: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
910: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
911: ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
912:
913: }
914: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
915: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
916: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
917: for (fc = 0; fc < Nc; ++fc) {
918: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
919: 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]);}
920: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
921: }
922: }
923: fieldOffset += Nb;
924: qc += Nc;
925: localDiff[field] += elemDiff;
926: }
927: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
928: }
929: DMRestoreLocalVector(dm, &localX);
930: MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
931: for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
932: PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
933: return(0);
934: }
936: /*@C
937: 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.
939: Input Parameters:
940: + dm - The DM
941: . time - The time
942: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
943: . ctxs - Optional array of contexts to pass to each function, or NULL.
944: - X - The coefficient vector u_h
946: Output Parameter:
947: . D - A Vec which holds the difference ||u - u_h||_2 for each cell
949: Level: developer
951: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
952: @*/
953: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
954: {
955: PetscSection section;
956: PetscQuadrature quad;
957: Vec localX;
958: PetscScalar *funcVal, *interpolant;
959: PetscReal *coords, *detJ, *J;
960: const PetscReal *quadPoints, *quadWeights;
961: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
962: PetscErrorCode ierr;
965: VecSet(D, 0.0);
966: DMGetDimension(dm, &dim);
967: DMGetCoordinateDim(dm, &coordDim);
968: DMGetSection(dm, §ion);
969: PetscSectionGetNumFields(section, &numFields);
970: DMGetLocalVector(dm, &localX);
971: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
972: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
973: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
974: for (field = 0; field < numFields; ++field) {
975: PetscObject obj;
976: PetscClassId id;
977: PetscInt Nc;
979: DMGetField(dm, field, &obj);
980: PetscObjectGetClassId(obj, &id);
981: if (id == PETSCFE_CLASSID) {
982: PetscFE fe = (PetscFE) obj;
984: PetscFEGetQuadrature(fe, &quad);
985: PetscFEGetNumComponents(fe, &Nc);
986: } else if (id == PETSCFV_CLASSID) {
987: PetscFV fv = (PetscFV) obj;
989: PetscFVGetQuadrature(fv, &quad);
990: PetscFVGetNumComponents(fv, &Nc);
991: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
992: numComponents += Nc;
993: }
994: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
995: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
996: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
997: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
998: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
999: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1000: for (c = cStart; c < cEnd; ++c) {
1001: PetscScalar *x = NULL;
1002: PetscScalar elemDiff = 0.0;
1003: PetscInt qc = 0;
1005: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
1006: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
1008: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1009: PetscObject obj;
1010: PetscClassId id;
1011: void * const ctx = ctxs ? ctxs[field] : NULL;
1012: PetscInt Nb, Nc, q, fc;
1014: DMGetField(dm, field, &obj);
1015: PetscObjectGetClassId(obj, &id);
1016: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1017: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1018: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1019: if (funcs[field]) {
1020: for (q = 0; q < Nq; ++q) {
1021: 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);
1022: (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
1023: if (ierr) {
1024: PetscErrorCode ierr2;
1025: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1026: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
1027: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1028:
1029: }
1030: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1031: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1032: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1033: for (fc = 0; fc < Nc; ++fc) {
1034: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1035: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
1036: }
1037: }
1038: }
1039: fieldOffset += Nb;
1040: qc += Nc;
1041: }
1042: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1043: VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1044: }
1045: PetscFree5(funcVal,interpolant,coords,detJ,J);
1046: DMRestoreLocalVector(dm, &localX);
1047: VecSqrtAbs(D);
1048: return(0);
1049: }
1051: /*@C
1052: DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec.
1054: Input Parameters:
1055: + dm - The DM
1056: - LocX - The coefficient vector u_h
1058: Output Parameter:
1059: . locC - A Vec which holds the Clement interpolant of the gradient
1061: Notes:
1062: Add citation to (Clement, 1975) and definition of the interpolant
1063: \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
1065: Level: developer
1067: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1068: @*/
1069: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1070: {
1071: DM_Plex *mesh = (DM_Plex *) dm->data;
1072: PetscInt debug = mesh->printFEM;
1073: DM dmC;
1074: PetscSection section;
1075: PetscQuadrature quad;
1076: PetscScalar *interpolant, *gradsum;
1077: PetscReal *coords, *detJ, *J, *invJ;
1078: const PetscReal *quadPoints, *quadWeights;
1079: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset;
1080: PetscErrorCode ierr;
1083: VecGetDM(locC, &dmC);
1084: VecSet(locC, 0.0);
1085: DMGetDimension(dm, &dim);
1086: DMGetCoordinateDim(dm, &coordDim);
1087: DMGetSection(dm, §ion);
1088: PetscSectionGetNumFields(section, &numFields);
1089: for (field = 0; field < numFields; ++field) {
1090: PetscObject obj;
1091: PetscClassId id;
1092: PetscInt Nc;
1094: DMGetField(dm, field, &obj);
1095: PetscObjectGetClassId(obj, &id);
1096: if (id == PETSCFE_CLASSID) {
1097: PetscFE fe = (PetscFE) obj;
1099: PetscFEGetQuadrature(fe, &quad);
1100: PetscFEGetNumComponents(fe, &Nc);
1101: } else if (id == PETSCFV_CLASSID) {
1102: PetscFV fv = (PetscFV) obj;
1104: PetscFVGetQuadrature(fv, &quad);
1105: PetscFVGetNumComponents(fv, &Nc);
1106: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1107: numComponents += Nc;
1108: }
1109: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1110: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1111: PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);
1112: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1113: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1114: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1115: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1116: for (v = vStart; v < vEnd; ++v) {
1117: PetscScalar volsum = 0.0;
1118: PetscInt *star = NULL;
1119: PetscInt starSize, st, d, fc;
1121: PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));
1122: DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1123: for (st = 0; st < starSize*2; st += 2) {
1124: const PetscInt cell = star[st];
1125: PetscScalar *grad = &gradsum[coordDim*numComponents];
1126: PetscScalar *x = NULL;
1127: PetscReal vol = 0.0;
1129: if ((cell < cStart) || (cell >= cEnd)) continue;
1130: DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);
1131: DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1132: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1133: PetscObject obj;
1134: PetscClassId id;
1135: PetscInt Nb, Nc, q, qc = 0;
1137: PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));
1138: DMGetField(dm, field, &obj);
1139: PetscObjectGetClassId(obj, &id);
1140: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1141: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1142: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1143: for (q = 0; q < Nq; ++q) {
1144: 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);
1145: if (ierr) {
1146: PetscErrorCode ierr2;
1147: ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1148: ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1149: ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2);
1150:
1151: }
1152: if (id == PETSCFE_CLASSID) {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);}
1153: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1154: for (fc = 0; fc < Nc; ++fc) {
1155: const PetscReal wt = quadWeights[q*qNc+qc+fc];
1157: for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q];
1158: }
1159: vol += quadWeights[q*qNc]*detJ[q];
1160: }
1161: fieldOffset += Nb;
1162: qc += Nc;
1163: }
1164: DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1165: for (fc = 0; fc < numComponents; ++fc) {
1166: for (d = 0; d < coordDim; ++d) {
1167: gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1168: }
1169: }
1170: volsum += vol;
1171: if (debug) {
1172: PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1173: for (fc = 0; fc < numComponents; ++fc) {
1174: for (d = 0; d < coordDim; ++d) {
1175: if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1176: PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1177: }
1178: }
1179: PetscPrintf(PETSC_COMM_SELF, "]\n");
1180: }
1181: }
1182: for (fc = 0; fc < numComponents; ++fc) {
1183: for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1184: }
1185: DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1186: DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1187: }
1188: PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);
1189: return(0);
1190: }
1192: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1193: {
1194: DM dmAux = NULL;
1195: PetscDS prob, probAux = NULL;
1196: PetscSection section, sectionAux;
1197: Vec locX, locA;
1198: PetscInt dim, numCells = cEnd - cStart, c, f;
1199: PetscBool useFVM = PETSC_FALSE;
1200: /* DS */
1201: PetscInt Nf, totDim, *uOff, *uOff_x, numConstants;
1202: PetscInt NfAux, totDimAux, *aOff;
1203: PetscScalar *u, *a;
1204: const PetscScalar *constants;
1205: /* Geometry */
1206: PetscFEGeom *cgeomFEM;
1207: DM dmGrad;
1208: PetscQuadrature affineQuad = NULL;
1209: Vec cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1210: PetscFVCellGeom *cgeomFVM;
1211: const PetscScalar *lgrad;
1212: PetscInt maxDegree;
1213: DMField coordField;
1214: IS cellIS;
1215: PetscErrorCode ierr;
1218: DMGetDS(dm, &prob);
1219: DMGetDimension(dm, &dim);
1220: DMGetSection(dm, §ion);
1221: PetscSectionGetNumFields(section, &Nf);
1222: /* Determine which discretizations we have */
1223: for (f = 0; f < Nf; ++f) {
1224: PetscObject obj;
1225: PetscClassId id;
1227: PetscDSGetDiscretization(prob, f, &obj);
1228: PetscObjectGetClassId(obj, &id);
1229: if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1230: }
1231: /* Get local solution with boundary values */
1232: DMGetLocalVector(dm, &locX);
1233: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1234: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1235: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1236: /* Read DS information */
1237: PetscDSGetTotalDimension(prob, &totDim);
1238: PetscDSGetComponentOffsets(prob, &uOff);
1239: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1240: ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1241: PetscDSGetConstants(prob, &numConstants, &constants);
1242: /* Read Auxiliary DS information */
1243: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1244: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1245: if (dmAux) {
1246: DMGetDS(dmAux, &probAux);
1247: PetscDSGetNumFields(probAux, &NfAux);
1248: DMGetSection(dmAux, §ionAux);
1249: PetscDSGetTotalDimension(probAux, &totDimAux);
1250: PetscDSGetComponentOffsets(probAux, &aOff);
1251: }
1252: /* Allocate data arrays */
1253: PetscCalloc1(numCells*totDim, &u);
1254: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1255: /* Read out geometry */
1256: DMGetCoordinateField(dm,&coordField);
1257: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1258: if (maxDegree <= 1) {
1259: DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1260: if (affineQuad) {
1261: DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1262: }
1263: }
1264: if (useFVM) {
1265: PetscFV fv = NULL;
1266: Vec grad;
1267: PetscInt fStart, fEnd;
1268: PetscBool compGrad;
1270: for (f = 0; f < Nf; ++f) {
1271: PetscObject obj;
1272: PetscClassId id;
1274: PetscDSGetDiscretization(prob, f, &obj);
1275: PetscObjectGetClassId(obj, &id);
1276: if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1277: }
1278: PetscFVGetComputeGradients(fv, &compGrad);
1279: PetscFVSetComputeGradients(fv, PETSC_TRUE);
1280: DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1281: DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1282: PetscFVSetComputeGradients(fv, compGrad);
1283: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1284: /* Reconstruct and limit cell gradients */
1285: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1286: DMGetGlobalVector(dmGrad, &grad);
1287: DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1288: /* Communicate gradient values */
1289: DMGetLocalVector(dmGrad, &locGrad);
1290: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1291: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1292: DMRestoreGlobalVector(dmGrad, &grad);
1293: /* Handle non-essential (e.g. outflow) boundary values */
1294: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1295: VecGetArrayRead(locGrad, &lgrad);
1296: }
1297: /* Read out data from inputs */
1298: for (c = cStart; c < cEnd; ++c) {
1299: PetscScalar *x = NULL;
1300: PetscInt i;
1302: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1303: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1304: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1305: if (dmAux) {
1306: DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1307: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1308: DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1309: }
1310: }
1311: /* Do integration for each field */
1312: for (f = 0; f < Nf; ++f) {
1313: PetscObject obj;
1314: PetscClassId id;
1315: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1317: PetscDSGetDiscretization(prob, f, &obj);
1318: PetscObjectGetClassId(obj, &id);
1319: if (id == PETSCFE_CLASSID) {
1320: PetscFE fe = (PetscFE) obj;
1321: PetscQuadrature q;
1322: PetscFEGeom *chunkGeom = NULL;
1323: PetscInt Nq, Nb;
1325: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1326: PetscFEGetQuadrature(fe, &q);
1327: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1328: PetscFEGetDimension(fe, &Nb);
1329: blockSize = Nb*Nq;
1330: batchSize = numBlocks * blockSize;
1331: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1332: numChunks = numCells / (numBatches*batchSize);
1333: Ne = numChunks*numBatches*batchSize;
1334: Nr = numCells % (numBatches*batchSize);
1335: offset = numCells - Nr;
1336: if (!affineQuad) {
1337: DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
1338: }
1339: PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1340: PetscFEIntegrate(fe, prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
1341: PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
1342: PetscFEIntegrate(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1343: PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
1344: if (!affineQuad) {
1345: PetscFEGeomDestroy(&cgeomFEM);
1346: }
1347: } else if (id == PETSCFV_CLASSID) {
1348: PetscInt foff;
1349: PetscPointFunc obj_func;
1350: PetscScalar lint;
1352: PetscDSGetObjective(prob, f, &obj_func);
1353: PetscDSGetFieldOffset(prob, f, &foff);
1354: if (obj_func) {
1355: for (c = 0; c < numCells; ++c) {
1356: PetscScalar *u_x;
1358: DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1359: 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);
1360: cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1361: }
1362: }
1363: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1364: }
1365: /* Cleanup data arrays */
1366: if (useFVM) {
1367: VecRestoreArrayRead(locGrad, &lgrad);
1368: VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1369: DMRestoreLocalVector(dmGrad, &locGrad);
1370: VecDestroy(&faceGeometryFVM);
1371: VecDestroy(&cellGeometryFVM);
1372: DMDestroy(&dmGrad);
1373: }
1374: if (dmAux) {PetscFree(a);}
1375: PetscFree(u);
1376: /* Cleanup */
1377: if (affineQuad) {
1378: PetscFEGeomDestroy(&cgeomFEM);
1379: }
1380: PetscQuadratureDestroy(&affineQuad);
1381: ISDestroy(&cellIS);
1382: DMRestoreLocalVector(dm, &locX);
1383: return(0);
1384: }
1386: /*@
1387: DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user
1389: Input Parameters:
1390: + dm - The mesh
1391: . X - Global input vector
1392: - user - The user context
1394: Output Parameter:
1395: . integral - Integral for each field
1397: Level: developer
1399: .seealso: DMPlexComputeResidualFEM()
1400: @*/
1401: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1402: {
1403: DM_Plex *mesh = (DM_Plex *) dm->data;
1404: PetscScalar *cintegral, *lintegral;
1405: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1412: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1413: DMGetNumFields(dm, &Nf);
1414: DMPlexGetVTKCellHeight(dm, &cellHeight);
1415: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1416: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1417: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1418: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1419: PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1420: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1421: /* Sum up values */
1422: for (cell = cStart; cell < cEnd; ++cell) {
1423: const PetscInt c = cell - cStart;
1425: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1426: for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1427: }
1428: MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1429: if (mesh->printFEM) {
1430: PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1431: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1432: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1433: }
1434: PetscFree2(lintegral, cintegral);
1435: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1436: return(0);
1437: }
1439: /*@
1440: DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user
1442: Input Parameters:
1443: + dm - The mesh
1444: . X - Global input vector
1445: - user - The user context
1447: Output Parameter:
1448: . integral - Cellwise integrals for each field
1450: Level: developer
1452: .seealso: DMPlexComputeResidualFEM()
1453: @*/
1454: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
1455: {
1456: DM_Plex *mesh = (DM_Plex *) dm->data;
1457: DM dmF;
1458: PetscSection sectionF;
1459: PetscScalar *cintegral, *af;
1460: PetscInt Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;
1467: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1468: DMGetNumFields(dm, &Nf);
1469: DMPlexGetVTKCellHeight(dm, &cellHeight);
1470: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1471: DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1472: cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1473: /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1474: PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
1475: DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1476: /* Put values in F*/
1477: VecGetDM(F, &dmF);
1478: DMGetSection(dmF, §ionF);
1479: VecGetArray(F, &af);
1480: for (cell = cStart; cell < cEnd; ++cell) {
1481: const PetscInt c = cell - cStart;
1482: PetscInt dof, off;
1484: if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1485: PetscSectionGetDof(sectionF, cell, &dof);
1486: PetscSectionGetOffset(sectionF, cell, &off);
1487: if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
1488: for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
1489: }
1490: VecRestoreArray(F, &af);
1491: PetscFree(cintegral);
1492: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1493: return(0);
1494: }
1496: static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS,
1497: void (*func)(PetscInt, PetscInt, PetscInt,
1498: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1499: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1500: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1501: PetscScalar *fintegral, void *user)
1502: {
1503: DM plex = NULL, plexA = NULL;
1504: PetscDS prob, probAux = NULL;
1505: PetscSection section, sectionAux = NULL;
1506: Vec locA = NULL;
1507: DMField coordField;
1508: PetscInt Nf, totDim, *uOff, *uOff_x;
1509: PetscInt NfAux = 0, totDimAux = 0, *aOff = NULL;
1510: PetscScalar *u, *a = NULL;
1511: const PetscScalar *constants;
1512: PetscInt numConstants, f;
1513: PetscErrorCode ierr;
1516: DMGetCoordinateField(dm, &coordField);
1517: DMConvert(dm, DMPLEX, &plex);
1518: DMGetDS(dm, &prob);
1519: DMGetDefaultSection(dm, §ion);
1520: PetscSectionGetNumFields(section, &Nf);
1521: /* Determine which discretizations we have */
1522: for (f = 0; f < Nf; ++f) {
1523: PetscObject obj;
1524: PetscClassId id;
1526: PetscDSGetDiscretization(prob, f, &obj);
1527: PetscObjectGetClassId(obj, &id);
1528: if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f);
1529: }
1530: /* Read DS information */
1531: PetscDSGetTotalDimension(prob, &totDim);
1532: PetscDSGetComponentOffsets(prob, &uOff);
1533: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1534: PetscDSGetConstants(prob, &numConstants, &constants);
1535: /* Read Auxiliary DS information */
1536: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1537: if (locA) {
1538: DM dmAux;
1540: VecGetDM(locA, &dmAux);
1541: DMConvert(dmAux, DMPLEX, &plexA);
1542: DMGetDS(dmAux, &probAux);
1543: PetscDSGetNumFields(probAux, &NfAux);
1544: DMGetDefaultSection(dmAux, §ionAux);
1545: PetscDSGetTotalDimension(probAux, &totDimAux);
1546: PetscDSGetComponentOffsets(probAux, &aOff);
1547: }
1548: /* Integrate over points */
1549: {
1550: PetscFEGeom *fgeom, *chunkGeom = NULL;
1551: PetscInt maxDegree;
1552: PetscQuadrature qGeom = NULL;
1553: const PetscInt *points;
1554: PetscInt numFaces, face, Nq, field;
1555: PetscInt numChunks, chunkSize, chunk, Nr, offset;
1557: ISGetLocalSize(pointIS, &numFaces);
1558: ISGetIndices(pointIS, &points);
1559: PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);
1560: DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
1561: for (field = 0; field < Nf; ++field) {
1562: PetscFE fe;
1564: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1565: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);}
1566: if (!qGeom) {
1567: PetscFEGetFaceQuadrature(fe, &qGeom);
1568: PetscObjectReference((PetscObject) qGeom);
1569: }
1570: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
1571: DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
1572: for (face = 0; face < numFaces; ++face) {
1573: const PetscInt point = points[face], *support, *cone;
1574: PetscScalar *x = NULL;
1575: PetscInt i, coneSize, faceLoc;
1577: DMPlexGetSupport(dm, point, &support);
1578: DMPlexGetConeSize(dm, support[0], &coneSize);
1579: DMPlexGetCone(dm, support[0], &cone);
1580: for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
1581: if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]);
1582: fgeom->face[face][0] = faceLoc;
1583: DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
1584: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
1585: DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
1586: if (locA) {
1587: PetscInt subp;
1588: DMPlexGetSubpoint(plexA, support[0], &subp);
1589: DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
1590: for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i];
1591: DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
1592: }
1593: }
1594: /* Get blocking */
1595: {
1596: PetscQuadrature q;
1597: PetscInt numBatches, batchSize, numBlocks, blockSize;
1598: PetscInt Nq, Nb;
1600: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1601: PetscFEGetQuadrature(fe, &q);
1602: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1603: PetscFEGetDimension(fe, &Nb);
1604: blockSize = Nb*Nq;
1605: batchSize = numBlocks * blockSize;
1606: chunkSize = numBatches*batchSize;
1607: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1608: numChunks = numFaces / chunkSize;
1609: Nr = numFaces % chunkSize;
1610: offset = numFaces - Nr;
1611: }
1612: /* Do integration for each field */
1613: for (chunk = 0; chunk < numChunks; ++chunk) {
1614: PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);
1615: PetscFEIntegrateBd(fe, prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);
1616: PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
1617: }
1618: PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);
1619: PetscFEIntegrateBd(fe, prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);
1620: PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);
1621: /* Cleanup data arrays */
1622: DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
1623: PetscQuadratureDestroy(&qGeom);
1624: PetscFree2(u, a);
1625: ISRestoreIndices(pointIS, &points);
1626: }
1627: }
1628: if (plex) {DMDestroy(&plex);}
1629: if (plexA) {DMDestroy(&plexA);}
1630: return(0);
1631: }
1633: /*@
1634: DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user
1636: Input Parameters:
1637: + dm - The mesh
1638: . X - Global input vector
1639: . label - The boundary DMLabel
1640: . numVals - The number of label values to use, or PETSC_DETERMINE for all values
1641: . vals - The label values to use, or PETSC_NULL for all values
1642: . func = The function to integrate along the boundary
1643: - user - The user context
1645: Output Parameter:
1646: . integral - Integral for each field
1648: Level: developer
1650: .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM()
1651: @*/
1652: PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[],
1653: void (*func)(PetscInt, PetscInt, PetscInt,
1654: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1655: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1656: PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1657: PetscScalar *integral, void *user)
1658: {
1659: Vec locX;
1660: PetscSection section;
1661: DMLabel depthLabel;
1662: IS facetIS;
1663: PetscInt dim, Nf, f, v;
1672: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1673: DMPlexGetDepthLabel(dm, &depthLabel);
1674: DMGetDimension(dm, &dim);
1675: DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
1676: DMGetDefaultSection(dm, §ion);
1677: PetscSectionGetNumFields(section, &Nf);
1678: /* Get local solution with boundary values */
1679: DMGetLocalVector(dm, &locX);
1680: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1681: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1682: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1683: /* Loop over label values */
1684: PetscMemzero(integral, Nf * sizeof(PetscScalar));
1685: for (v = 0; v < numVals; ++v) {
1686: IS pointIS;
1687: PetscInt numFaces, face;
1688: PetscScalar *fintegral;
1690: DMLabelGetStratumIS(label, vals[v], &pointIS);
1691: if (!pointIS) continue; /* No points with that id on this process */
1692: {
1693: IS isectIS;
1695: /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
1696: ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);
1697: ISDestroy(&pointIS);
1698: pointIS = isectIS;
1699: }
1700: ISGetLocalSize(pointIS, &numFaces);
1701: PetscCalloc1(numFaces*Nf, &fintegral);
1702: DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);
1703: /* Sum point contributions into integral */
1704: for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f];
1705: PetscFree(fintegral);
1706: ISDestroy(&pointIS);
1707: }
1708: DMRestoreLocalVector(dm, &locX);
1709: ISDestroy(&facetIS);
1710: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1711: return(0);
1712: }
1714: /*@
1715: DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1717: Input Parameters:
1718: + dmf - The fine mesh
1719: . dmc - The coarse mesh
1720: - user - The user context
1722: Output Parameter:
1723: . In - The interpolation matrix
1725: Level: developer
1727: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1728: @*/
1729: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1730: {
1731: DM_Plex *mesh = (DM_Plex *) dmc->data;
1732: const char *name = "Interpolator";
1733: PetscDS prob;
1734: PetscFE *feRef;
1735: PetscFV *fvRef;
1736: PetscSection fsection, fglobalSection;
1737: PetscSection csection, cglobalSection;
1738: PetscScalar *elemMat;
1739: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1740: PetscInt cTotDim, rTotDim = 0;
1741: PetscErrorCode ierr;
1744: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1745: DMGetDimension(dmf, &dim);
1746: DMGetSection(dmf, &fsection);
1747: DMGetGlobalSection(dmf, &fglobalSection);
1748: DMGetSection(dmc, &csection);
1749: DMGetGlobalSection(dmc, &cglobalSection);
1750: PetscSectionGetNumFields(fsection, &Nf);
1751: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1752: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1753: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1754: DMGetDS(dmf, &prob);
1755: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1756: for (f = 0; f < Nf; ++f) {
1757: PetscObject obj;
1758: PetscClassId id;
1759: PetscInt rNb = 0, Nc = 0;
1761: PetscDSGetDiscretization(prob, f, &obj);
1762: PetscObjectGetClassId(obj, &id);
1763: if (id == PETSCFE_CLASSID) {
1764: PetscFE fe = (PetscFE) obj;
1766: PetscFERefine(fe, &feRef[f]);
1767: PetscFEGetDimension(feRef[f], &rNb);
1768: PetscFEGetNumComponents(fe, &Nc);
1769: } else if (id == PETSCFV_CLASSID) {
1770: PetscFV fv = (PetscFV) obj;
1771: PetscDualSpace Q;
1773: PetscFVRefine(fv, &fvRef[f]);
1774: PetscFVGetDualSpace(fvRef[f], &Q);
1775: PetscDualSpaceGetDimension(Q, &rNb);
1776: PetscFVGetNumComponents(fv, &Nc);
1777: }
1778: rTotDim += rNb;
1779: }
1780: PetscDSGetTotalDimension(prob, &cTotDim);
1781: PetscMalloc1(rTotDim*cTotDim,&elemMat);
1782: PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1783: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1784: PetscDualSpace Qref;
1785: PetscQuadrature f;
1786: const PetscReal *qpoints, *qweights;
1787: PetscReal *points;
1788: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
1790: /* Compose points from all dual basis functionals */
1791: if (feRef[fieldI]) {
1792: PetscFEGetDualSpace(feRef[fieldI], &Qref);
1793: PetscFEGetNumComponents(feRef[fieldI], &Nc);
1794: } else {
1795: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1796: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1797: }
1798: PetscDualSpaceGetDimension(Qref, &fpdim);
1799: for (i = 0; i < fpdim; ++i) {
1800: PetscDualSpaceGetFunctional(Qref, i, &f);
1801: PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1802: npoints += Np;
1803: }
1804: PetscMalloc1(npoints*dim,&points);
1805: for (i = 0, k = 0; i < fpdim; ++i) {
1806: PetscDualSpaceGetFunctional(Qref, i, &f);
1807: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1808: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1809: }
1811: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1812: PetscObject obj;
1813: PetscClassId id;
1814: PetscReal *B;
1815: PetscInt NcJ = 0, cpdim = 0, j, qNc;
1817: PetscDSGetDiscretization(prob, fieldJ, &obj);
1818: PetscObjectGetClassId(obj, &id);
1819: if (id == PETSCFE_CLASSID) {
1820: PetscFE fe = (PetscFE) obj;
1822: /* Evaluate basis at points */
1823: PetscFEGetNumComponents(fe, &NcJ);
1824: PetscFEGetDimension(fe, &cpdim);
1825: /* For now, fields only interpolate themselves */
1826: if (fieldI == fieldJ) {
1827: 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);
1828: PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1829: for (i = 0, k = 0; i < fpdim; ++i) {
1830: PetscDualSpaceGetFunctional(Qref, i, &f);
1831: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1832: 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);
1833: for (p = 0; p < Np; ++p, ++k) {
1834: for (j = 0; j < cpdim; ++j) {
1835: /*
1836: cTotDim: Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
1837: offsetI, offsetJ: Offsets into the larger element interpolation matrix for different fields
1838: fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
1839: qNC, Nc, Ncj, c: Number of components in this field
1840: Np, p: Number of quad points in the fine grid functional i
1841: k: i*Np + p, overall point number for the interpolation
1842: */
1843: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1844: }
1845: }
1846: }
1847: PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1848: }
1849: } else if (id == PETSCFV_CLASSID) {
1850: PetscFV fv = (PetscFV) obj;
1852: /* Evaluate constant function at points */
1853: PetscFVGetNumComponents(fv, &NcJ);
1854: cpdim = 1;
1855: /* For now, fields only interpolate themselves */
1856: if (fieldI == fieldJ) {
1857: 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);
1858: for (i = 0, k = 0; i < fpdim; ++i) {
1859: PetscDualSpaceGetFunctional(Qref, i, &f);
1860: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1861: 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);
1862: for (p = 0; p < Np; ++p, ++k) {
1863: for (j = 0; j < cpdim; ++j) {
1864: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c];
1865: }
1866: }
1867: }
1868: }
1869: }
1870: offsetJ += cpdim;
1871: }
1872: offsetI += fpdim;
1873: PetscFree(points);
1874: }
1875: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1876: /* Preallocate matrix */
1877: {
1878: Mat preallocator;
1879: PetscScalar *vals;
1880: PetscInt *cellCIndices, *cellFIndices;
1881: PetscInt locRows, locCols, cell;
1883: MatGetLocalSize(In, &locRows, &locCols);
1884: MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1885: MatSetType(preallocator, MATPREALLOCATOR);
1886: MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1887: MatSetUp(preallocator);
1888: PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1889: for (cell = cStart; cell < cEnd; ++cell) {
1890: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1891: MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1892: }
1893: PetscFree3(vals,cellCIndices,cellFIndices);
1894: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1895: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1896: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1897: MatDestroy(&preallocator);
1898: }
1899: /* Fill matrix */
1900: MatZeroEntries(In);
1901: for (c = cStart; c < cEnd; ++c) {
1902: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1903: }
1904: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1905: PetscFree2(feRef,fvRef);
1906: PetscFree(elemMat);
1907: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1908: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1909: if (mesh->printFEM) {
1910: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1911: MatChop(In, 1.0e-10);
1912: MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1913: }
1914: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1915: return(0);
1916: }
1918: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
1919: {
1920: SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
1921: }
1923: /*@
1924: DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.
1926: Input Parameters:
1927: + dmf - The fine mesh
1928: . dmc - The coarse mesh
1929: - user - The user context
1931: Output Parameter:
1932: . In - The interpolation matrix
1934: Level: developer
1936: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1937: @*/
1938: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1939: {
1940: DM_Plex *mesh = (DM_Plex *) dmf->data;
1941: const char *name = "Interpolator";
1942: PetscDS prob;
1943: PetscSection fsection, csection, globalFSection, globalCSection;
1944: PetscHSetIJ ht;
1945: PetscLayout rLayout;
1946: PetscInt *dnz, *onz;
1947: PetscInt locRows, rStart, rEnd;
1948: PetscReal *x, *v0, *J, *invJ, detJ;
1949: PetscReal *v0c, *Jc, *invJc, detJc;
1950: PetscScalar *elemMat;
1951: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
1955: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1956: DMGetCoordinateDim(dmc, &dim);
1957: DMGetDS(dmc, &prob);
1958: PetscDSGetRefCoordArrays(prob, &x, NULL);
1959: PetscDSGetNumFields(prob, &Nf);
1960: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1961: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1962: DMGetSection(dmf, &fsection);
1963: DMGetGlobalSection(dmf, &globalFSection);
1964: DMGetSection(dmc, &csection);
1965: DMGetGlobalSection(dmc, &globalCSection);
1966: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1967: PetscDSGetTotalDimension(prob, &totDim);
1968: PetscMalloc1(totDim, &elemMat);
1970: MatGetLocalSize(In, &locRows, NULL);
1971: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1972: PetscLayoutSetLocalSize(rLayout, locRows);
1973: PetscLayoutSetBlockSize(rLayout, 1);
1974: PetscLayoutSetUp(rLayout);
1975: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1976: PetscLayoutDestroy(&rLayout);
1977: PetscCalloc2(locRows,&dnz,locRows,&onz);
1978: PetscHSetIJCreate(&ht);
1979: for (field = 0; field < Nf; ++field) {
1980: PetscObject obj;
1981: PetscClassId id;
1982: PetscDualSpace Q = NULL;
1983: PetscQuadrature f;
1984: const PetscReal *qpoints;
1985: PetscInt Nc, Np, fpdim, i, d;
1987: PetscDSGetDiscretization(prob, field, &obj);
1988: PetscObjectGetClassId(obj, &id);
1989: if (id == PETSCFE_CLASSID) {
1990: PetscFE fe = (PetscFE) obj;
1992: PetscFEGetDualSpace(fe, &Q);
1993: PetscFEGetNumComponents(fe, &Nc);
1994: } else if (id == PETSCFV_CLASSID) {
1995: PetscFV fv = (PetscFV) obj;
1997: PetscFVGetDualSpace(fv, &Q);
1998: Nc = 1;
1999: }
2000: PetscDualSpaceGetDimension(Q, &fpdim);
2001: /* For each fine grid cell */
2002: for (cell = cStart; cell < cEnd; ++cell) {
2003: PetscInt *findices, *cindices;
2004: PetscInt numFIndices, numCIndices;
2006: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2007: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2008: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2009: for (i = 0; i < fpdim; ++i) {
2010: Vec pointVec;
2011: PetscScalar *pV;
2012: PetscSF coarseCellSF = NULL;
2013: const PetscSFNode *coarseCells;
2014: PetscInt numCoarseCells, q, c;
2016: /* Get points from the dual basis functional quadrature */
2017: PetscDualSpaceGetFunctional(Q, i, &f);
2018: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2019: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2020: VecSetBlockSize(pointVec, dim);
2021: VecGetArray(pointVec, &pV);
2022: for (q = 0; q < Np; ++q) {
2023: const PetscReal xi0[3] = {-1., -1., -1.};
2025: /* Transform point to real space */
2026: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2027: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2028: }
2029: VecRestoreArray(pointVec, &pV);
2030: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2031: /* OPT: Pack all quad points from fine cell */
2032: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2033: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2034: /* Update preallocation info */
2035: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2036: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2037: {
2038: PetscHashIJKey key;
2039: PetscBool missing;
2041: key.i = findices[i];
2042: if (key.i >= 0) {
2043: /* Get indices for coarse elements */
2044: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2045: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2046: for (c = 0; c < numCIndices; ++c) {
2047: key.j = cindices[c];
2048: if (key.j < 0) continue;
2049: PetscHSetIJQueryAdd(ht, key, &missing);
2050: if (missing) {
2051: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2052: else ++onz[key.i-rStart];
2053: }
2054: }
2055: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2056: }
2057: }
2058: }
2059: PetscSFDestroy(&coarseCellSF);
2060: VecDestroy(&pointVec);
2061: }
2062: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2063: }
2064: }
2065: PetscHSetIJDestroy(&ht);
2066: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
2067: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2068: PetscFree2(dnz,onz);
2069: for (field = 0; field < Nf; ++field) {
2070: PetscObject obj;
2071: PetscClassId id;
2072: PetscDualSpace Q = NULL;
2073: PetscQuadrature f;
2074: const PetscReal *qpoints, *qweights;
2075: PetscInt Nc, qNc, Np, fpdim, i, d;
2077: PetscDSGetDiscretization(prob, field, &obj);
2078: PetscObjectGetClassId(obj, &id);
2079: if (id == PETSCFE_CLASSID) {
2080: PetscFE fe = (PetscFE) obj;
2082: PetscFEGetDualSpace(fe, &Q);
2083: PetscFEGetNumComponents(fe, &Nc);
2084: } else if (id == PETSCFV_CLASSID) {
2085: PetscFV fv = (PetscFV) obj;
2087: PetscFVGetDualSpace(fv, &Q);
2088: Nc = 1;
2089: } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field);
2090: PetscDualSpaceGetDimension(Q, &fpdim);
2091: /* For each fine grid cell */
2092: for (cell = cStart; cell < cEnd; ++cell) {
2093: PetscInt *findices, *cindices;
2094: PetscInt numFIndices, numCIndices;
2096: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2097: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2098: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2099: for (i = 0; i < fpdim; ++i) {
2100: Vec pointVec;
2101: PetscScalar *pV;
2102: PetscSF coarseCellSF = NULL;
2103: const PetscSFNode *coarseCells;
2104: PetscInt numCoarseCells, cpdim, q, c, j;
2106: /* Get points from the dual basis functional quadrature */
2107: PetscDualSpaceGetFunctional(Q, i, &f);
2108: PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
2109: 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);
2110: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2111: VecSetBlockSize(pointVec, dim);
2112: VecGetArray(pointVec, &pV);
2113: for (q = 0; q < Np; ++q) {
2114: const PetscReal xi0[3] = {-1., -1., -1.};
2116: /* Transform point to real space */
2117: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2118: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2119: }
2120: VecRestoreArray(pointVec, &pV);
2121: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2122: /* OPT: Read this out from preallocation information */
2123: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2124: /* Update preallocation info */
2125: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2126: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2127: VecGetArray(pointVec, &pV);
2128: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2129: PetscReal pVReal[3];
2130: const PetscReal xi0[3] = {-1., -1., -1.};
2132: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2133: /* Transform points from real space to coarse reference space */
2134: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2135: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2136: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
2138: if (id == PETSCFE_CLASSID) {
2139: PetscFE fe = (PetscFE) obj;
2140: PetscReal *B;
2142: /* Evaluate coarse basis on contained point */
2143: PetscFEGetDimension(fe, &cpdim);
2144: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2145: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2146: /* Get elemMat entries by multiplying by weight */
2147: for (j = 0; j < cpdim; ++j) {
2148: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
2149: }
2150: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2151: } else {
2152: cpdim = 1;
2153: for (j = 0; j < cpdim; ++j) {
2154: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
2155: }
2156: }
2157: /* Update interpolator */
2158: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2159: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2160: MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
2161: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2162: }
2163: VecRestoreArray(pointVec, &pV);
2164: PetscSFDestroy(&coarseCellSF);
2165: VecDestroy(&pointVec);
2166: }
2167: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2168: }
2169: }
2170: PetscFree3(v0,J,invJ);
2171: PetscFree3(v0c,Jc,invJc);
2172: PetscFree(elemMat);
2173: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2174: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2175: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2176: return(0);
2177: }
2179: /*@
2180: DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM.
2182: Input Parameters:
2183: + dmf - The fine mesh
2184: . dmc - The coarse mesh
2185: - user - The user context
2187: Output Parameter:
2188: . mass - The mass matrix
2190: Level: developer
2192: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2193: @*/
2194: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
2195: {
2196: DM_Plex *mesh = (DM_Plex *) dmf->data;
2197: const char *name = "Mass Matrix";
2198: PetscDS prob;
2199: PetscSection fsection, csection, globalFSection, globalCSection;
2200: PetscHSetIJ ht;
2201: PetscLayout rLayout;
2202: PetscInt *dnz, *onz;
2203: PetscInt locRows, rStart, rEnd;
2204: PetscReal *x, *v0, *J, *invJ, detJ;
2205: PetscReal *v0c, *Jc, *invJc, detJc;
2206: PetscScalar *elemMat;
2207: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
2211: DMGetCoordinateDim(dmc, &dim);
2212: DMGetDS(dmc, &prob);
2213: PetscDSGetRefCoordArrays(prob, &x, NULL);
2214: PetscDSGetNumFields(prob, &Nf);
2215: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2216: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2217: DMGetSection(dmf, &fsection);
2218: DMGetGlobalSection(dmf, &globalFSection);
2219: DMGetSection(dmc, &csection);
2220: DMGetGlobalSection(dmc, &globalCSection);
2221: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2222: PetscDSGetTotalDimension(prob, &totDim);
2223: PetscMalloc1(totDim, &elemMat);
2225: MatGetLocalSize(mass, &locRows, NULL);
2226: PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
2227: PetscLayoutSetLocalSize(rLayout, locRows);
2228: PetscLayoutSetBlockSize(rLayout, 1);
2229: PetscLayoutSetUp(rLayout);
2230: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2231: PetscLayoutDestroy(&rLayout);
2232: PetscCalloc2(locRows,&dnz,locRows,&onz);
2233: PetscHSetIJCreate(&ht);
2234: for (field = 0; field < Nf; ++field) {
2235: PetscObject obj;
2236: PetscClassId id;
2237: PetscQuadrature quad;
2238: const PetscReal *qpoints;
2239: PetscInt Nq, Nc, i, d;
2241: PetscDSGetDiscretization(prob, field, &obj);
2242: PetscObjectGetClassId(obj, &id);
2243: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
2244: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2245: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
2246: /* For each fine grid cell */
2247: for (cell = cStart; cell < cEnd; ++cell) {
2248: Vec pointVec;
2249: PetscScalar *pV;
2250: PetscSF coarseCellSF = NULL;
2251: const PetscSFNode *coarseCells;
2252: PetscInt numCoarseCells, q, c;
2253: PetscInt *findices, *cindices;
2254: PetscInt numFIndices, numCIndices;
2256: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2257: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2258: /* Get points from the quadrature */
2259: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2260: VecSetBlockSize(pointVec, dim);
2261: VecGetArray(pointVec, &pV);
2262: for (q = 0; q < Nq; ++q) {
2263: const PetscReal xi0[3] = {-1., -1., -1.};
2265: /* Transform point to real space */
2266: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2267: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2268: }
2269: VecRestoreArray(pointVec, &pV);
2270: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2271: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2272: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2273: /* Update preallocation info */
2274: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2275: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2276: {
2277: PetscHashIJKey key;
2278: PetscBool missing;
2280: for (i = 0; i < numFIndices; ++i) {
2281: key.i = findices[i];
2282: if (key.i >= 0) {
2283: /* Get indices for coarse elements */
2284: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2285: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2286: for (c = 0; c < numCIndices; ++c) {
2287: key.j = cindices[c];
2288: if (key.j < 0) continue;
2289: PetscHSetIJQueryAdd(ht, key, &missing);
2290: if (missing) {
2291: if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2292: else ++onz[key.i-rStart];
2293: }
2294: }
2295: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2296: }
2297: }
2298: }
2299: }
2300: PetscSFDestroy(&coarseCellSF);
2301: VecDestroy(&pointVec);
2302: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2303: }
2304: }
2305: PetscHSetIJDestroy(&ht);
2306: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2307: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2308: PetscFree2(dnz,onz);
2309: for (field = 0; field < Nf; ++field) {
2310: PetscObject obj;
2311: PetscClassId id;
2312: PetscQuadrature quad;
2313: PetscReal *Bfine;
2314: const PetscReal *qpoints, *qweights;
2315: PetscInt Nq, Nc, i, d;
2317: PetscDSGetDiscretization(prob, field, &obj);
2318: PetscObjectGetClassId(obj, &id);
2319: if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);}
2320: else {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2321: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
2322: /* For each fine grid cell */
2323: for (cell = cStart; cell < cEnd; ++cell) {
2324: Vec pointVec;
2325: PetscScalar *pV;
2326: PetscSF coarseCellSF = NULL;
2327: const PetscSFNode *coarseCells;
2328: PetscInt numCoarseCells, cpdim, q, c, j;
2329: PetscInt *findices, *cindices;
2330: PetscInt numFIndices, numCIndices;
2332: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2333: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2334: /* Get points from the quadrature */
2335: VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2336: VecSetBlockSize(pointVec, dim);
2337: VecGetArray(pointVec, &pV);
2338: for (q = 0; q < Nq; ++q) {
2339: const PetscReal xi0[3] = {-1., -1., -1.};
2341: /* Transform point to real space */
2342: CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2343: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2344: }
2345: VecRestoreArray(pointVec, &pV);
2346: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2347: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2348: /* Update matrix */
2349: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2350: if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2351: VecGetArray(pointVec, &pV);
2352: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2353: PetscReal pVReal[3];
2354: const PetscReal xi0[3] = {-1., -1., -1.};
2357: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2358: /* Transform points from real space to coarse reference space */
2359: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2360: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2361: CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);
2363: if (id == PETSCFE_CLASSID) {
2364: PetscFE fe = (PetscFE) obj;
2365: PetscReal *B;
2367: /* Evaluate coarse basis on contained point */
2368: PetscFEGetDimension(fe, &cpdim);
2369: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2370: /* Get elemMat entries by multiplying by weight */
2371: for (i = 0; i < numFIndices; ++i) {
2372: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2373: for (j = 0; j < cpdim; ++j) {
2374: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2375: }
2376: /* Update interpolator */
2377: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2378: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2379: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2380: }
2381: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2382: } else {
2383: cpdim = 1;
2384: for (i = 0; i < numFIndices; ++i) {
2385: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2386: for (j = 0; j < cpdim; ++j) {
2387: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2388: }
2389: /* Update interpolator */
2390: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2391: PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);
2392: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2393: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2394: }
2395: }
2396: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2397: }
2398: VecRestoreArray(pointVec, &pV);
2399: PetscSFDestroy(&coarseCellSF);
2400: VecDestroy(&pointVec);
2401: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2402: }
2403: }
2404: PetscFree3(v0,J,invJ);
2405: PetscFree3(v0c,Jc,invJc);
2406: PetscFree(elemMat);
2407: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2408: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2409: return(0);
2410: }
2412: /*@
2413: DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns
2415: Input Parameters:
2416: + dmc - The coarse mesh
2417: - dmf - The fine mesh
2418: - user - The user context
2420: Output Parameter:
2421: . sc - The mapping
2423: Level: developer
2425: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2426: @*/
2427: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
2428: {
2429: PetscDS prob;
2430: PetscFE *feRef;
2431: PetscFV *fvRef;
2432: Vec fv, cv;
2433: IS fis, cis;
2434: PetscSection fsection, fglobalSection, csection, cglobalSection;
2435: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
2436: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
2437: PetscBool *needAvg;
2441: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2442: DMGetDimension(dmf, &dim);
2443: DMGetSection(dmf, &fsection);
2444: DMGetGlobalSection(dmf, &fglobalSection);
2445: DMGetSection(dmc, &csection);
2446: DMGetGlobalSection(dmc, &cglobalSection);
2447: PetscSectionGetNumFields(fsection, &Nf);
2448: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2449: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2450: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2451: DMGetDS(dmc, &prob);
2452: PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);
2453: for (f = 0; f < Nf; ++f) {
2454: PetscObject obj;
2455: PetscClassId id;
2456: PetscInt fNb = 0, Nc = 0;
2458: PetscDSGetDiscretization(prob, f, &obj);
2459: PetscObjectGetClassId(obj, &id);
2460: if (id == PETSCFE_CLASSID) {
2461: PetscFE fe = (PetscFE) obj;
2462: PetscSpace sp;
2463: PetscInt maxDegree;
2465: PetscFERefine(fe, &feRef[f]);
2466: PetscFEGetDimension(feRef[f], &fNb);
2467: PetscFEGetNumComponents(fe, &Nc);
2468: PetscFEGetBasisSpace(fe, &sp);
2469: PetscSpaceGetDegree(sp, NULL, &maxDegree);
2470: if (!maxDegree) needAvg[f] = PETSC_TRUE;
2471: } else if (id == PETSCFV_CLASSID) {
2472: PetscFV fv = (PetscFV) obj;
2473: PetscDualSpace Q;
2475: PetscFVRefine(fv, &fvRef[f]);
2476: PetscFVGetDualSpace(fvRef[f], &Q);
2477: PetscDualSpaceGetDimension(Q, &fNb);
2478: PetscFVGetNumComponents(fv, &Nc);
2479: needAvg[f] = PETSC_TRUE;
2480: }
2481: fTotDim += fNb;
2482: }
2483: PetscDSGetTotalDimension(prob, &cTotDim);
2484: PetscMalloc1(cTotDim,&cmap);
2485: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
2486: PetscFE feC;
2487: PetscFV fvC;
2488: PetscDualSpace QF, QC;
2489: PetscInt order = -1, NcF, NcC, fpdim, cpdim;
2491: if (feRef[field]) {
2492: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
2493: PetscFEGetNumComponents(feC, &NcC);
2494: PetscFEGetNumComponents(feRef[field], &NcF);
2495: PetscFEGetDualSpace(feRef[field], &QF);
2496: PetscDualSpaceGetOrder(QF, &order);
2497: PetscDualSpaceGetDimension(QF, &fpdim);
2498: PetscFEGetDualSpace(feC, &QC);
2499: PetscDualSpaceGetDimension(QC, &cpdim);
2500: } else {
2501: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
2502: PetscFVGetNumComponents(fvC, &NcC);
2503: PetscFVGetNumComponents(fvRef[field], &NcF);
2504: PetscFVGetDualSpace(fvRef[field], &QF);
2505: PetscDualSpaceGetDimension(QF, &fpdim);
2506: PetscFVGetDualSpace(fvC, &QC);
2507: PetscDualSpaceGetDimension(QC, &cpdim);
2508: }
2509: 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);
2510: for (c = 0; c < cpdim; ++c) {
2511: PetscQuadrature cfunc;
2512: const PetscReal *cqpoints, *cqweights;
2513: PetscInt NqcC, NpC;
2514: PetscBool found = PETSC_FALSE;
2516: PetscDualSpaceGetFunctional(QC, c, &cfunc);
2517: PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
2518: if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC);
2519: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
2520: for (f = 0; f < fpdim; ++f) {
2521: PetscQuadrature ffunc;
2522: const PetscReal *fqpoints, *fqweights;
2523: PetscReal sum = 0.0;
2524: PetscInt NqcF, NpF;
2526: PetscDualSpaceGetFunctional(QF, f, &ffunc);
2527: PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
2528: if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF);
2529: if (NpC != NpF) continue;
2530: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
2531: if (sum > 1.0e-9) continue;
2532: for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
2533: if (sum < 1.0e-9) continue;
2534: cmap[offsetC+c] = offsetF+f;
2535: found = PETSC_TRUE;
2536: break;
2537: }
2538: if (!found) {
2539: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
2540: if (fvRef[field] || (feRef[field] && order == 0)) {
2541: cmap[offsetC+c] = offsetF+0;
2542: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
2543: }
2544: }
2545: offsetC += cpdim;
2546: offsetF += fpdim;
2547: }
2548: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
2549: PetscFree3(feRef,fvRef,needAvg);
2551: DMGetGlobalVector(dmf, &fv);
2552: DMGetGlobalVector(dmc, &cv);
2553: VecGetOwnershipRange(cv, &startC, &endC);
2554: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
2555: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
2556: PetscMalloc1(m,&cindices);
2557: PetscMalloc1(m,&findices);
2558: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
2559: for (c = cStart; c < cEnd; ++c) {
2560: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
2561: for (d = 0; d < cTotDim; ++d) {
2562: if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
2563: 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]]);
2564: cindices[cellCIndices[d]-startC] = cellCIndices[d];
2565: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
2566: }
2567: }
2568: PetscFree(cmap);
2569: PetscFree2(cellCIndices,cellFIndices);
2571: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
2572: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
2573: VecScatterCreate(cv, cis, fv, fis, sc);
2574: ISDestroy(&cis);
2575: ISDestroy(&fis);
2576: DMRestoreGlobalVector(dmf, &fv);
2577: DMRestoreGlobalVector(dmc, &cv);
2578: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2579: return(0);
2580: }
2582: PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
2583: {
2584: char composeStr[33] = {0};
2585: PetscObjectId id;
2586: PetscContainer container;
2587: PetscErrorCode ierr;
2590: PetscObjectGetId((PetscObject)quad,&id);
2591: PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);
2592: PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
2593: if (container) {
2594: PetscContainerGetPointer(container, (void **) geom);
2595: } else {
2596: DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
2597: PetscContainerCreate(PETSC_COMM_SELF,&container);
2598: PetscContainerSetPointer(container, (void *) *geom);
2599: PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
2600: PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
2601: PetscContainerDestroy(&container);
2602: }
2603: return(0);
2604: }
2606: PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
2607: {
2609: *geom = NULL;
2610: return(0);
2611: }
2613: /*
2614: 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
2616: X - The local solution vector
2617: X_t - The local solution time derviative vector, or NULL
2618: */
2619: PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS,
2620: PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx)
2621: {
2622: DM_Plex *mesh = (DM_Plex *) dm->data;
2623: const char *name = "Jacobian", *nameP = "JacobianPre";
2624: DM dmAux = NULL;
2625: PetscDS prob, probAux = NULL;
2626: PetscSection sectionAux = NULL;
2627: Vec A;
2628: DMField coordField;
2629: PetscFEGeom *cgeomFEM;
2630: PetscQuadrature qGeom = NULL;
2631: Mat J = Jac, JP = JacP;
2632: PetscScalar *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL;
2633: PetscBool hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE;
2634: const PetscInt *cells;
2635: PetscInt Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0;
2636: PetscErrorCode ierr;
2639: CHKMEMQ;
2640: ISGetLocalSize(cellIS, &numCells);
2641: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
2642: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
2643: DMGetDS(dm, &prob);
2644: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2645: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2646: if (dmAux) {
2647: DMGetDefaultSection(dmAux, §ionAux);
2648: DMGetDS(dmAux, &probAux);
2649: }
2650: /* Get flags */
2651: PetscDSGetNumFields(prob, &Nf);
2652: DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);
2653: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2654: PetscObject disc;
2655: PetscClassId id;
2656: PetscDSGetDiscretization(prob, fieldI, &disc);
2657: PetscObjectGetClassId(disc, &id);
2658: if (id == PETSCFE_CLASSID) {isFE[fieldI] = PETSC_TRUE;}
2659: else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;}
2660: }
2661: PetscDSHasJacobian(prob, &hasJac);
2662: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
2663: PetscDSHasDynamicJacobian(prob, &hasDyn);
2664: assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE;
2665: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
2666: PetscObjectTypeCompare((PetscObject) Jac, MATIS, &isMatIS);
2667: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
2668: /* Setup input data and temp arrays (should be DMGetWorkArray) */
2669: if (isMatISP || isMatISP) {DMPlexGetSubdomainSection(dm, &globalSection);}
2670: if (isMatIS) {MatISGetLocalMat(Jac, &J);}
2671: if (isMatISP) {MatISGetLocalMat(JacP, &JP);}
2672: if (hasFV) {MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);} /* No allocated space for FV stuff, so ignore the zero entries */
2673: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2674: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2675: PetscDSGetTotalDimension(prob, &totDim);
2676: if (probAux) {PetscDSGetTotalDimension(probAux, &totDimAux);}
2677: CHKMEMQ;
2678: /* Compute batch sizes */
2679: if (isFE[0]) {
2680: PetscFE fe;
2681: PetscQuadrature q;
2682: PetscInt numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb;
2684: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
2685: PetscFEGetQuadrature(fe, &q);
2686: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
2687: PetscFEGetDimension(fe, &Nb);
2688: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2689: blockSize = Nb*numQuadPoints;
2690: batchSize = numBlocks * blockSize;
2691: chunkSize = numBatches * batchSize;
2692: numChunks = numCells / chunkSize + numCells % chunkSize;
2693: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2694: } else {
2695: chunkSize = numCells;
2696: numChunks = 1;
2697: }
2698: /* Get work space */
2699: wsz = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize;
2700: DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);
2701: PetscMemzero(work, wsz * sizeof(PetscScalar));
2702: off = 0;
2703: u = X ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
2704: u_t = X_t ? (sz = chunkSize*totDim, off += sz, work+off-sz) : NULL;
2705: a = dmAux ? (sz = chunkSize*totDimAux, off += sz, work+off-sz) : NULL;
2706: elemMat = hasJac ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
2707: elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
2708: elemMatD = hasDyn ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
2709: if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz);
2710: /* Setup geometry */
2711: DMGetCoordinateField(dm, &coordField);
2712: DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
2713: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);}
2714: if (!qGeom) {
2715: PetscFE fe;
2717: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
2718: PetscFEGetQuadrature(fe, &qGeom);
2719: PetscObjectReference((PetscObject) qGeom);
2720: }
2721: DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
2722: /* Compute volume integrals */
2723: if (assembleJac) {MatZeroEntries(J);}
2724: MatZeroEntries(JP);
2725: for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) {
2726: const PetscInt Ncell = PetscMin(chunkSize, numCells - offCell);
2727: PetscInt c;
2729: /* Extract values */
2730: for (c = 0; c < Ncell; ++c) {
2731: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
2732: PetscScalar *x = NULL, *x_t = NULL;
2733: PetscInt i;
2735: if (X) {
2736: DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
2737: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
2738: DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
2739: }
2740: if (X_t) {
2741: DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
2742: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
2743: DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
2744: }
2745: if (dmAux) {
2746: DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);
2747: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
2748: DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);
2749: }
2750: }
2751: CHKMEMQ;
2752: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2753: PetscFE fe;
2754: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2755: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2756: if (hasJac) {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);}
2757: if (hasPrec) {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);}
2758: if (hasDyn) {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);}
2759: }
2760: /* For finite volume, add the identity */
2761: if (!isFE[fieldI]) {
2762: PetscFV fv;
2763: PetscInt eOffset = 0, Nc, fc, foff;
2765: PetscDSGetFieldOffset(prob, fieldI, &foff);
2766: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
2767: PetscFVGetNumComponents(fv, &Nc);
2768: for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) {
2769: for (fc = 0; fc < Nc; ++fc) {
2770: const PetscInt i = foff + fc;
2771: if (hasJac) {elemMat [eOffset+i*totDim+i] = 1.0;}
2772: if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;}
2773: }
2774: }
2775: }
2776: }
2777: CHKMEMQ;
2778: /* Add contribution from X_t */
2779: if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
2780: /* Insert values into matrix */
2781: for (c = 0; c < Ncell; ++c) {
2782: const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
2783: if (mesh->printFEM > 1) {
2784: if (hasJac) {DMPrintCellMatrix(cell, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2785: if (hasPrec) {DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
2786: }
2787: if (assembleJac) {DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);}
2788: DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2789: }
2790: CHKMEMQ;
2791: }
2792: /* Cleanup */
2793: DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
2794: PetscQuadratureDestroy(&qGeom);
2795: if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
2796: DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);
2797: 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);
2798: /* Compute boundary integrals */
2799: /* DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx); */
2800: /* Assemble matrix */
2801: if (assembleJac) {MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);}
2802: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
2803: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2804: CHKMEMQ;
2805: return(0);
2806: }