Actual source code: plexfem.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petsc/private/petscfvimpl.h>
7: /*@
8: DMPlexGetScale - Get the scale for the specified fundamental unit
10: Not collective
12: Input Arguments:
13: + dm - the DM
14: - unit - The SI unit
16: Output Argument:
17: . scale - The value used to scale all quantities with this unit
19: Level: advanced
21: .seealso: DMPlexSetScale(), PetscUnit
22: @*/
23: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
24: {
25: DM_Plex *mesh = (DM_Plex*) dm->data;
30: *scale = mesh->scale[unit];
31: return(0);
32: }
34: /*@
35: DMPlexSetScale - Set the scale for the specified fundamental unit
37: Not collective
39: Input Arguments:
40: + dm - the DM
41: . unit - The SI unit
42: - scale - The value used to scale all quantities with this unit
44: Level: advanced
46: .seealso: DMPlexGetScale(), PetscUnit
47: @*/
48: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
49: {
50: DM_Plex *mesh = (DM_Plex*) dm->data;
54: mesh->scale[unit] = scale;
55: return(0);
56: }
58: PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
59: {
60: switch (i) {
61: case 0:
62: switch (j) {
63: case 0: return 0;
64: case 1:
65: switch (k) {
66: case 0: return 0;
67: case 1: return 0;
68: case 2: return 1;
69: }
70: case 2:
71: switch (k) {
72: case 0: return 0;
73: case 1: return -1;
74: case 2: return 0;
75: }
76: }
77: case 1:
78: switch (j) {
79: case 0:
80: switch (k) {
81: case 0: return 0;
82: case 1: return 0;
83: case 2: return -1;
84: }
85: case 1: return 0;
86: case 2:
87: switch (k) {
88: case 0: return 1;
89: case 1: return 0;
90: case 2: return 0;
91: }
92: }
93: case 2:
94: switch (j) {
95: case 0:
96: switch (k) {
97: case 0: return 0;
98: case 1: return 1;
99: case 2: return 0;
100: }
101: case 1:
102: switch (k) {
103: case 0: return -1;
104: case 1: return 0;
105: case 2: return 0;
106: }
107: case 2: return 0;
108: }
109: }
110: return 0;
111: }
113: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
114: {
115: PetscInt *ctxInt = (PetscInt *) ctx;
116: PetscInt dim2 = ctxInt[0];
117: PetscInt d = ctxInt[1];
118: PetscInt i, j, k = dim > 2 ? d - dim : d;
121: if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
122: for (i = 0; i < dim; i++) mode[i] = 0.;
123: if (d < dim) {
124: mode[d] = 1.;
125: } else {
126: for (i = 0; i < dim; i++) {
127: for (j = 0; j < dim; j++) {
128: mode[j] += epsilon(i, j, k)*X[i];
129: }
130: }
131: }
132: return(0);
133: }
135: /*@C
136: DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates
138: Collective on DM
140: Input Arguments:
141: . dm - the DM
143: Output Argument:
144: . sp - the null space
146: Note: This is necessary to take account of Dirichlet conditions on the displacements
148: Level: advanced
150: .seealso: MatNullSpaceCreate()
151: @*/
152: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
153: {
154: MPI_Comm comm;
155: Vec mode[6];
156: PetscSection section, globalSection;
157: PetscInt dim, dimEmbed, n, m, d, i, j;
161: PetscObjectGetComm((PetscObject)dm,&comm);
162: DMGetDimension(dm, &dim);
163: DMGetCoordinateDim(dm, &dimEmbed);
164: if (dim == 1) {
165: MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
166: return(0);
167: }
168: DMGetDefaultSection(dm, §ion);
169: DMGetDefaultGlobalSection(dm, &globalSection);
170: PetscSectionGetConstrainedStorageSize(globalSection, &n);
171: m = (dim*(dim+1))/2;
172: VecCreate(comm, &mode[0]);
173: VecSetSizes(mode[0], n, PETSC_DETERMINE);
174: VecSetUp(mode[0]);
175: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
176: for (d = 0; d < m; d++) {
177: PetscInt ctx[2];
178: PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
179: void *voidctx = (void *) (&ctx[0]);
181: ctx[0] = dimEmbed;
182: ctx[1] = d;
183: DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
184: }
185: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
186: /* Orthonormalize system */
187: for (i = dim; i < m; ++i) {
188: PetscScalar dots[6];
190: VecMDot(mode[i], i, mode, dots);
191: for (j = 0; j < i; ++j) dots[j] *= -1.0;
192: VecMAXPY(mode[i], i, dots, mode);
193: VecNormalize(mode[i], NULL);
194: }
195: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
196: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
197: return(0);
198: }
200: /*@
201: DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
202: are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
203: evaluating the dual space basis of that point. A basis function is associated with the point in its
204: transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
205: projection height, which is set with this function. By default, the maximum projection height is zero, which means
206: that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior
207: basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
209: Input Parameters:
210: + dm - the DMPlex object
211: - height - the maximum projection height >= 0
213: Level: advanced
215: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
216: @*/
217: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
218: {
219: DM_Plex *plex = (DM_Plex *) dm->data;
223: plex->maxProjectionHeight = height;
224: return(0);
225: }
227: /*@
228: DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
229: DMPlexProjectXXXLocal() functions.
231: Input Parameters:
232: . dm - the DMPlex object
234: Output Parameters:
235: . height - the maximum projection height
237: Level: intermediate
239: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
240: @*/
241: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
242: {
243: DM_Plex *plex = (DM_Plex *) dm->data;
247: *height = plex->maxProjectionHeight;
248: return(0);
249: }
251: /*@C
252: DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector
254: Input Parameters:
255: + dm - The DM, with a PetscDS that matches the problem being constrained
256: . time - The time
257: . field - The field to constrain
258: . Nc - The number of constrained field components, or 0 for all components
259: . comps - An array of constrained component numbers, or NULL for all components
260: . label - The DMLabel defining constrained points
261: . numids - The number of DMLabel ids for constrained points
262: . ids - An array of ids for constrained points
263: . func - A pointwise function giving boundary values
264: - ctx - An optional user context for bcFunc
266: Output Parameter:
267: . locX - A local vector to receives the boundary values
269: Level: developer
271: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
272: @*/
273: 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)
274: {
275: PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
276: void **ctxs;
277: PetscInt numFields;
278: PetscErrorCode ierr;
281: DMGetNumFields(dm, &numFields);
282: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
283: funcs[field] = func;
284: ctxs[field] = ctx;
285: DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
286: PetscFree2(funcs,ctxs);
287: return(0);
288: }
290: /*@C
291: DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector
293: Input Parameters:
294: + dm - The DM, with a PetscDS that matches the problem being constrained
295: . time - The time
296: . locU - A local vector with the input solution values
297: . field - The field to constrain
298: . Nc - The number of constrained field components, or 0 for all components
299: . comps - An array of constrained component numbers, or NULL for all components
300: . label - The DMLabel defining constrained points
301: . numids - The number of DMLabel ids for constrained points
302: . ids - An array of ids for constrained points
303: . func - A pointwise function giving boundary values
304: - ctx - An optional user context for bcFunc
306: Output Parameter:
307: . locX - A local vector to receives the boundary values
309: Level: developer
311: .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary()
312: @*/
313: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
314: void (*func)(PetscInt, PetscInt, PetscInt,
315: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
316: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
317: PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
318: PetscScalar[]),
319: void *ctx, Vec locX)
320: {
321: void (**funcs)(PetscInt, PetscInt, PetscInt,
322: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
323: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
324: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
325: void **ctxs;
326: PetscInt numFields;
327: PetscErrorCode ierr;
330: DMGetNumFields(dm, &numFields);
331: PetscCalloc2(numFields,&funcs,numFields,&ctxs);
332: funcs[field] = func;
333: ctxs[field] = ctx;
334: DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
335: PetscFree2(funcs,ctxs);
336: return(0);
337: }
339: /*@C
340: DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector
342: Input Parameters:
343: + dm - The DM, with a PetscDS that matches the problem being constrained
344: . time - The time
345: . faceGeometry - A vector with the FVM face geometry information
346: . cellGeometry - A vector with the FVM cell geometry information
347: . Grad - A vector with the FVM cell gradient information
348: . field - The field to constrain
349: . Nc - The number of constrained field components, or 0 for all components
350: . comps - An array of constrained component numbers, or NULL for all components
351: . label - The DMLabel defining constrained points
352: . numids - The number of DMLabel ids for constrained points
353: . ids - An array of ids for constrained points
354: . func - A pointwise function giving boundary values
355: - ctx - An optional user context for bcFunc
357: Output Parameter:
358: . locX - A local vector to receives the boundary values
360: Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary()
362: Level: developer
364: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
365: @*/
366: 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[],
367: PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
368: {
369: PetscDS prob;
370: PetscSF sf;
371: DM dmFace, dmCell, dmGrad;
372: const PetscScalar *facegeom, *cellgeom = NULL, *grad;
373: const PetscInt *leaves;
374: PetscScalar *x, *fx;
375: PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i;
376: PetscErrorCode ierr, ierru = 0;
379: DMGetPointSF(dm, &sf);
380: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
381: nleaves = PetscMax(0, nleaves);
382: DMGetDimension(dm, &dim);
383: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
384: DMGetDS(dm, &prob);
385: VecGetDM(faceGeometry, &dmFace);
386: VecGetArrayRead(faceGeometry, &facegeom);
387: if (cellGeometry) {
388: VecGetDM(cellGeometry, &dmCell);
389: VecGetArrayRead(cellGeometry, &cellgeom);
390: }
391: if (Grad) {
392: PetscFV fv;
394: PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
395: VecGetDM(Grad, &dmGrad);
396: VecGetArrayRead(Grad, &grad);
397: PetscFVGetNumComponents(fv, &pdim);
398: DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
399: }
400: VecGetArray(locX, &x);
401: for (i = 0; i < numids; ++i) {
402: IS faceIS;
403: const PetscInt *faces;
404: PetscInt numFaces, f;
406: DMLabelGetStratumIS(label, ids[i], &faceIS);
407: if (!faceIS) continue; /* No points with that id on this process */
408: ISGetLocalSize(faceIS, &numFaces);
409: ISGetIndices(faceIS, &faces);
410: for (f = 0; f < numFaces; ++f) {
411: const PetscInt face = faces[f], *cells;
412: PetscFVFaceGeom *fg;
414: if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
415: PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
416: if (loc >= 0) continue;
417: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
418: DMPlexGetSupport(dm, face, &cells);
419: if (Grad) {
420: PetscFVCellGeom *cg;
421: PetscScalar *cx, *cgrad;
422: PetscScalar *xG;
423: PetscReal dx[3];
424: PetscInt d;
426: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
427: DMPlexPointLocalRead(dm, cells[0], x, &cx);
428: DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
429: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
430: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
431: for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
432: ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
433: if (ierru) {
434: ISRestoreIndices(faceIS, &faces);
435: ISDestroy(&faceIS);
436: goto cleanup;
437: }
438: } else {
439: PetscScalar *xI;
440: PetscScalar *xG;
442: DMPlexPointLocalRead(dm, cells[0], x, &xI);
443: DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
444: ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
445: if (ierru) {
446: ISRestoreIndices(faceIS, &faces);
447: ISDestroy(&faceIS);
448: goto cleanup;
449: }
450: }
451: }
452: ISRestoreIndices(faceIS, &faces);
453: ISDestroy(&faceIS);
454: }
455: cleanup:
456: VecRestoreArray(locX, &x);
457: if (Grad) {
458: DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
459: VecRestoreArrayRead(Grad, &grad);
460: }
461: if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
462: VecRestoreArrayRead(faceGeometry, &facegeom);
463: CHKERRQ(ierru);
464: return(0);
465: }
467: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
468: {
469: PetscInt numBd, b;
473: PetscDSGetNumBoundary(dm->prob, &numBd);
474: for (b = 0; b < numBd; ++b) {
475: DMBoundaryConditionType type;
476: const char *labelname;
477: DMLabel label;
478: PetscInt field, Nc;
479: const PetscInt *comps;
480: PetscObject obj;
481: PetscClassId id;
482: void (*func)(void);
483: PetscInt numids;
484: const PetscInt *ids;
485: void *ctx;
487: DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
488: if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
489: DMGetLabel(dm, labelname, &label);
490: DMGetField(dm, field, &obj);
491: PetscObjectGetClassId(obj, &id);
492: if (id == PETSCFE_CLASSID) {
493: switch (type) {
494: /* for FEM, there is no insertion to be done for non-essential boundary conditions */
495: case DM_BC_ESSENTIAL:
496: DMPlexLabelAddCells(dm,label);
497: DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
498: DMPlexLabelClearCells(dm,label);
499: break;
500: case DM_BC_ESSENTIAL_FIELD:
501: DMPlexLabelAddCells(dm,label);
502: DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
503: (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
504: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
505: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
506: DMPlexLabelClearCells(dm,label);
507: break;
508: default: break;
509: }
510: } else if (id == PETSCFV_CLASSID) {
511: if (!faceGeomFVM) continue;
512: DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
513: (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
514: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
515: }
516: return(0);
517: }
519: /*@
520: DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector
522: Input Parameters:
523: + dm - The DM
524: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
525: . time - The time
526: . faceGeomFVM - Face geometry data for FV discretizations
527: . cellGeomFVM - Cell geometry data for FV discretizations
528: - gradFVM - Gradient reconstruction data for FV discretizations
530: Output Parameters:
531: . locX - Solution updated with boundary values
533: Level: developer
535: .seealso: DMProjectFunctionLabelLocal()
536: @*/
537: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
538: {
547: PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
548: return(0);
549: }
551: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
552: {
553: const PetscInt debug = 0;
554: PetscSection section;
555: PetscQuadrature quad;
556: Vec localX;
557: PetscScalar *funcVal, *interpolant;
558: PetscReal *coords, *detJ, *J;
559: PetscReal localDiff = 0.0;
560: const PetscReal *quadWeights;
561: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
562: PetscErrorCode ierr;
565: DMGetDimension(dm, &dim);
566: DMGetCoordinateDim(dm, &coordDim);
567: DMGetDefaultSection(dm, §ion);
568: PetscSectionGetNumFields(section, &numFields);
569: DMGetLocalVector(dm, &localX);
570: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
571: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
572: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
573: for (field = 0; field < numFields; ++field) {
574: PetscObject obj;
575: PetscClassId id;
576: PetscInt Nc;
578: DMGetField(dm, field, &obj);
579: PetscObjectGetClassId(obj, &id);
580: if (id == PETSCFE_CLASSID) {
581: PetscFE fe = (PetscFE) obj;
583: PetscFEGetQuadrature(fe, &quad);
584: PetscFEGetNumComponents(fe, &Nc);
585: } else if (id == PETSCFV_CLASSID) {
586: PetscFV fv = (PetscFV) obj;
588: PetscFVGetQuadrature(fv, &quad);
589: PetscFVGetNumComponents(fv, &Nc);
590: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
591: numComponents += Nc;
592: }
593: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
594: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
595: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
596: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
597: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
598: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
599: for (c = cStart; c < cEnd; ++c) {
600: PetscScalar *x = NULL;
601: PetscReal elemDiff = 0.0;
602: PetscInt qc = 0;
604: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
605: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
607: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
608: PetscObject obj;
609: PetscClassId id;
610: void * const ctx = ctxs ? ctxs[field] : NULL;
611: PetscInt Nb, Nc, q, fc;
613: DMGetField(dm, field, &obj);
614: PetscObjectGetClassId(obj, &id);
615: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
616: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
617: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
618: if (debug) {
619: char title[1024];
620: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
621: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
622: }
623: for (q = 0; q < Nq; ++q) {
624: 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);
625: (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
626: if (ierr) {
627: PetscErrorCode ierr2;
628: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
629: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
630: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
631:
632: }
633: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
634: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
635: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
636: for (fc = 0; fc < Nc; ++fc) {
637: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
638: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
639: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
640: }
641: }
642: fieldOffset += Nb;
643: qc += Nc;
644: }
645: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
646: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
647: localDiff += elemDiff;
648: }
649: PetscFree5(funcVal,interpolant,coords,detJ,J);
650: DMRestoreLocalVector(dm, &localX);
651: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
652: *diff = PetscSqrtReal(*diff);
653: return(0);
654: }
656: 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)
657: {
658: const PetscInt debug = 0;
659: PetscSection section;
660: PetscQuadrature quad;
661: Vec localX;
662: PetscScalar *funcVal, *interpolant;
663: const PetscReal *quadPoints, *quadWeights;
664: PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ;
665: PetscReal localDiff = 0.0;
666: PetscInt dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
667: PetscErrorCode ierr;
670: DMGetDimension(dm, &dim);
671: DMGetCoordinateDim(dm, &coordDim);
672: DMGetDefaultSection(dm, §ion);
673: PetscSectionGetNumFields(section, &numFields);
674: DMGetLocalVector(dm, &localX);
675: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
676: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
677: for (field = 0; field < numFields; ++field) {
678: PetscFE fe;
679: PetscInt Nc;
681: DMGetField(dm, field, (PetscObject *) &fe);
682: PetscFEGetQuadrature(fe, &quad);
683: PetscFEGetNumComponents(fe, &Nc);
684: numComponents += Nc;
685: }
686: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
687: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
688: /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
689: PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
690: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
691: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
692: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
693: for (c = cStart; c < cEnd; ++c) {
694: PetscScalar *x = NULL;
695: PetscReal elemDiff = 0.0;
696: PetscInt qc = 0;
698: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);
699: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
701: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
702: PetscFE fe;
703: void * const ctx = ctxs ? ctxs[field] : NULL;
704: PetscReal *basisDer;
705: PetscInt Nb, Nc, q, fc;
707: DMGetField(dm, field, (PetscObject *) &fe);
708: PetscFEGetDimension(fe, &Nb);
709: PetscFEGetNumComponents(fe, &Nc);
710: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
711: if (debug) {
712: char title[1024];
713: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
714: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
715: }
716: for (q = 0; q < Nq; ++q) {
717: 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);
718: (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
719: if (ierr) {
720: PetscErrorCode ierr2;
721: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
722: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
723: ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
724:
725: }
726: PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
727: for (fc = 0; fc < Nc; ++fc) {
728: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
729: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
730: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
731: }
732: }
733: fieldOffset += Nb;
734: qc += Nc;
735: }
736: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
737: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
738: localDiff += elemDiff;
739: }
740: PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
741: DMRestoreLocalVector(dm, &localX);
742: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
743: *diff = PetscSqrtReal(*diff);
744: return(0);
745: }
747: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
748: {
749: const PetscInt debug = 0;
750: PetscSection section;
751: PetscQuadrature quad;
752: Vec localX;
753: PetscScalar *funcVal, *interpolant;
754: PetscReal *coords, *detJ, *J;
755: PetscReal *localDiff;
756: const PetscReal *quadPoints, *quadWeights;
757: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
758: PetscErrorCode ierr;
761: DMGetDimension(dm, &dim);
762: DMGetCoordinateDim(dm, &coordDim);
763: DMGetDefaultSection(dm, §ion);
764: PetscSectionGetNumFields(section, &numFields);
765: DMGetLocalVector(dm, &localX);
766: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
767: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
768: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
769: for (field = 0; field < numFields; ++field) {
770: PetscObject obj;
771: PetscClassId id;
772: PetscInt Nc;
774: DMGetField(dm, field, &obj);
775: PetscObjectGetClassId(obj, &id);
776: if (id == PETSCFE_CLASSID) {
777: PetscFE fe = (PetscFE) obj;
779: PetscFEGetQuadrature(fe, &quad);
780: PetscFEGetNumComponents(fe, &Nc);
781: } else if (id == PETSCFV_CLASSID) {
782: PetscFV fv = (PetscFV) obj;
784: PetscFVGetQuadrature(fv, &quad);
785: PetscFVGetNumComponents(fv, &Nc);
786: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
787: numComponents += Nc;
788: }
789: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
790: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
791: PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
792: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
793: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
794: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
795: for (c = cStart; c < cEnd; ++c) {
796: PetscScalar *x = NULL;
797: PetscInt qc = 0;
799: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
800: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
802: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
803: PetscObject obj;
804: PetscClassId id;
805: void * const ctx = ctxs ? ctxs[field] : NULL;
806: PetscInt Nb, Nc, q, fc;
808: PetscReal elemDiff = 0.0;
810: DMGetField(dm, field, &obj);
811: PetscObjectGetClassId(obj, &id);
812: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
813: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
814: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
815: if (debug) {
816: char title[1024];
817: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
818: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
819: }
820: for (q = 0; q < Nq; ++q) {
821: 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);
822: (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
823: if (ierr) {
824: PetscErrorCode ierr2;
825: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
826: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
827: ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
828:
829: }
830: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
831: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
832: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
833: for (fc = 0; fc < Nc; ++fc) {
834: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
835: 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]);}
836: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
837: }
838: }
839: fieldOffset += Nb;
840: qc += Nc;
841: localDiff[field] += elemDiff;
842: }
843: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
844: }
845: DMRestoreLocalVector(dm, &localX);
846: MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
847: for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
848: PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
849: return(0);
850: }
852: /*@C
853: 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.
855: Input Parameters:
856: + dm - The DM
857: . time - The time
858: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
859: . ctxs - Optional array of contexts to pass to each function, or NULL.
860: - X - The coefficient vector u_h
862: Output Parameter:
863: . D - A Vec which holds the difference ||u - u_h||_2 for each cell
865: Level: developer
867: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
868: @*/
869: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
870: {
871: PetscSection section;
872: PetscQuadrature quad;
873: Vec localX;
874: PetscScalar *funcVal, *interpolant;
875: PetscReal *coords, *detJ, *J;
876: const PetscReal *quadPoints, *quadWeights;
877: PetscInt dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
878: PetscErrorCode ierr;
881: VecSet(D, 0.0);
882: DMGetDimension(dm, &dim);
883: DMGetCoordinateDim(dm, &coordDim);
884: DMGetDefaultSection(dm, §ion);
885: PetscSectionGetNumFields(section, &numFields);
886: DMGetLocalVector(dm, &localX);
887: DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
888: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
889: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
890: for (field = 0; field < numFields; ++field) {
891: PetscObject obj;
892: PetscClassId id;
893: PetscInt Nc;
895: DMGetField(dm, field, &obj);
896: PetscObjectGetClassId(obj, &id);
897: if (id == PETSCFE_CLASSID) {
898: PetscFE fe = (PetscFE) obj;
900: PetscFEGetQuadrature(fe, &quad);
901: PetscFEGetNumComponents(fe, &Nc);
902: } else if (id == PETSCFV_CLASSID) {
903: PetscFV fv = (PetscFV) obj;
905: PetscFVGetQuadrature(fv, &quad);
906: PetscFVGetNumComponents(fv, &Nc);
907: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
908: numComponents += Nc;
909: }
910: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
911: if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
912: PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
913: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
914: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
915: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
916: for (c = cStart; c < cEnd; ++c) {
917: PetscScalar *x = NULL;
918: PetscScalar elemDiff = 0.0;
919: PetscInt qc = 0;
921: DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
922: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
924: for (field = 0, fieldOffset = 0; field < numFields; ++field) {
925: PetscObject obj;
926: PetscClassId id;
927: void * const ctx = ctxs ? ctxs[field] : NULL;
928: PetscInt Nb, Nc, q, fc;
930: DMGetField(dm, field, &obj);
931: PetscObjectGetClassId(obj, &id);
932: if (id == PETSCFE_CLASSID) {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
933: else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
934: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
935: if (funcs[field]) {
936: for (q = 0; q < Nq; ++q) {
937: 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);
938: (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx);
939: if (ierr) {
940: PetscErrorCode ierr2;
941: ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
942: ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
943: ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
944:
945: }
946: if (id == PETSCFE_CLASSID) {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
947: else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
948: else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
949: for (fc = 0; fc < Nc; ++fc) {
950: const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
951: elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
952: }
953: }
954: }
955: fieldOffset += Nb;
956: qc += Nc;
957: }
958: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
959: VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
960: }
961: PetscFree5(funcVal,interpolant,coords,detJ,J);
962: DMRestoreLocalVector(dm, &localX);
963: VecSqrtAbs(D);
964: return(0);
965: }
967: /*@
968: DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user
970: Input Parameters:
971: + dm - The mesh
972: . X - Local input vector
973: - user - The user context
975: Output Parameter:
976: . integral - Local integral for each field
978: Level: developer
980: .seealso: DMPlexComputeResidualFEM()
981: @*/
982: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
983: {
984: DM_Plex *mesh = (DM_Plex *) dm->data;
985: DM dmAux, dmGrad;
986: Vec localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
987: PetscDS prob, probAux = NULL;
988: PetscSection section, sectionAux;
989: PetscFV fvm = NULL;
990: PetscFECellGeom *cgeomFEM;
991: PetscFVFaceGeom *fgeomFVM;
992: PetscFVCellGeom *cgeomFVM;
993: PetscScalar *u, *a = NULL;
994: const PetscScalar *constants, *lgrad;
995: PetscReal *lintegral;
996: PetscInt *uOff, *uOff_x, *aOff = NULL;
997: PetscInt dim, numConstants, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c;
998: PetscInt totDim, totDimAux;
999: PetscBool useFVM = PETSC_FALSE;
1000: PetscErrorCode ierr;
1003: PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1004: DMGetLocalVector(dm, &localX);
1005: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);
1006: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1007: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1008: DMGetDimension(dm, &dim);
1009: DMGetDefaultSection(dm, §ion);
1010: DMGetDS(dm, &prob);
1011: PetscDSGetTotalDimension(prob, &totDim);
1012: PetscDSGetComponentOffsets(prob, &uOff);
1013: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1014: PetscDSGetConstants(prob, &numConstants, &constants);
1015: PetscSectionGetNumFields(section, &Nf);
1016: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1017: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1018: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1019: numCells = cEnd - cStart;
1020: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1021: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1022: if (dmAux) {
1023: DMGetDS(dmAux, &probAux);
1024: PetscDSGetNumFields(probAux, &NfAux);
1025: DMGetDefaultSection(dmAux, §ionAux);
1026: PetscDSGetTotalDimension(probAux, &totDimAux);
1027: PetscDSGetComponentOffsets(probAux, &aOff);
1028: PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);
1029: }
1030: for (f = 0; f < Nf; ++f) {
1031: PetscObject obj;
1032: PetscClassId id;
1034: PetscDSGetDiscretization(prob, f, &obj);
1035: PetscObjectGetClassId(obj, &id);
1036: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1037: }
1038: if (useFVM) {
1039: Vec grad;
1040: PetscInt fStart, fEnd;
1041: PetscBool compGrad;
1043: PetscFVGetComputeGradients(fvm, &compGrad);
1044: PetscFVSetComputeGradients(fvm, PETSC_TRUE);
1045: DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1046: DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1047: PetscFVSetComputeGradients(fvm, compGrad);
1048: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1049: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1050: /* Reconstruct and limit cell gradients */
1051: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1052: DMGetGlobalVector(dmGrad, &grad);
1053: DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);
1054: /* Communicate gradient values */
1055: DMGetLocalVector(dmGrad, &locGrad);
1056: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1057: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1058: DMRestoreGlobalVector(dmGrad, &grad);
1059: /* Handle non-essential (e.g. outflow) boundary values */
1060: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1061: VecGetArrayRead(locGrad, &lgrad);
1062: }
1063: PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);
1064: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1065: for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1066: for (c = cStart; c < cEnd; ++c) {
1067: PetscScalar *x = NULL;
1068: PetscInt i;
1070: DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);
1071: if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c);
1072: DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1073: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1074: DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1075: if (dmAux) {
1076: DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1077: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1078: DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1079: }
1080: }
1081: for (f = 0; f < Nf; ++f) {
1082: PetscObject obj;
1083: PetscClassId id;
1084: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1086: PetscDSGetDiscretization(prob, f, &obj);
1087: PetscObjectGetClassId(obj, &id);
1088: if (id == PETSCFE_CLASSID) {
1089: PetscFE fe = (PetscFE) obj;
1090: PetscQuadrature q;
1091: PetscInt Nq, Nb;
1093: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1094: PetscFEGetQuadrature(fe, &q);
1095: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1096: PetscFEGetDimension(fe, &Nb);
1097: blockSize = Nb*Nq;
1098: batchSize = numBlocks * blockSize;
1099: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1100: numChunks = numCells / (numBatches*batchSize);
1101: Ne = numChunks*numBatches*batchSize;
1102: Nr = numCells % (numBatches*batchSize);
1103: offset = numCells - Nr;
1104: PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);
1105: PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1106: } else if (id == PETSCFV_CLASSID) {
1107: /* PetscFV fv = (PetscFV) obj; */
1108: PetscInt foff;
1109: PetscPointFunc obj_func;
1110: PetscScalar lint;
1112: PetscDSGetObjective(prob, f, &obj_func);
1113: PetscDSGetFieldOffset(prob, f, &foff);
1114: if (obj_func) {
1115: for (c = 0; c < numCells; ++c) {
1116: PetscScalar *u_x;
1118: DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1119: obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, a, NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint);
1120: lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1121: }
1122: }
1123: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1124: }
1125: if (useFVM) {
1126: VecRestoreArrayRead(locGrad, &lgrad);
1127: VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1128: VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1129: DMRestoreLocalVector(dmGrad, &locGrad);
1130: VecDestroy(&faceGeometryFVM);
1131: VecDestroy(&cellGeometryFVM);
1132: DMDestroy(&dmGrad);
1133: }
1134: if (dmAux) {PetscFree(a);}
1135: if (mesh->printFEM) {
1136: PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1137: for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1138: PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1139: }
1140: DMRestoreLocalVector(dm, &localX);
1141: MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1142: PetscFree3(lintegral,u,cgeomFEM);
1143: PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1144: return(0);
1145: }
1147: /*@
1148: DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1150: Input Parameters:
1151: + dmf - The fine mesh
1152: . dmc - The coarse mesh
1153: - user - The user context
1155: Output Parameter:
1156: . In - The interpolation matrix
1158: Level: developer
1160: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1161: @*/
1162: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1163: {
1164: DM_Plex *mesh = (DM_Plex *) dmc->data;
1165: const char *name = "Interpolator";
1166: PetscDS prob;
1167: PetscFE *feRef;
1168: PetscFV *fvRef;
1169: PetscSection fsection, fglobalSection;
1170: PetscSection csection, cglobalSection;
1171: PetscScalar *elemMat;
1172: PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1173: PetscInt cTotDim, rTotDim = 0;
1174: PetscErrorCode ierr;
1177: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1178: DMGetDimension(dmf, &dim);
1179: DMGetDefaultSection(dmf, &fsection);
1180: DMGetDefaultGlobalSection(dmf, &fglobalSection);
1181: DMGetDefaultSection(dmc, &csection);
1182: DMGetDefaultGlobalSection(dmc, &cglobalSection);
1183: PetscSectionGetNumFields(fsection, &Nf);
1184: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1185: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1186: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1187: DMGetDS(dmf, &prob);
1188: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1189: for (f = 0; f < Nf; ++f) {
1190: PetscObject obj;
1191: PetscClassId id;
1192: PetscInt rNb = 0, Nc = 0;
1194: PetscDSGetDiscretization(prob, f, &obj);
1195: PetscObjectGetClassId(obj, &id);
1196: if (id == PETSCFE_CLASSID) {
1197: PetscFE fe = (PetscFE) obj;
1199: PetscFERefine(fe, &feRef[f]);
1200: PetscFEGetDimension(feRef[f], &rNb);
1201: PetscFEGetNumComponents(fe, &Nc);
1202: } else if (id == PETSCFV_CLASSID) {
1203: PetscFV fv = (PetscFV) obj;
1204: PetscDualSpace Q;
1206: PetscFVRefine(fv, &fvRef[f]);
1207: PetscFVGetDualSpace(fvRef[f], &Q);
1208: PetscDualSpaceGetDimension(Q, &rNb);
1209: PetscFVGetNumComponents(fv, &Nc);
1210: }
1211: rTotDim += rNb;
1212: }
1213: PetscDSGetTotalDimension(prob, &cTotDim);
1214: PetscMalloc1(rTotDim*cTotDim,&elemMat);
1215: PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1216: for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1217: PetscDualSpace Qref;
1218: PetscQuadrature f;
1219: const PetscReal *qpoints, *qweights;
1220: PetscReal *points;
1221: PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d;
1223: /* Compose points from all dual basis functionals */
1224: if (feRef[fieldI]) {
1225: PetscFEGetDualSpace(feRef[fieldI], &Qref);
1226: PetscFEGetNumComponents(feRef[fieldI], &Nc);
1227: } else {
1228: PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1229: PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1230: }
1231: PetscDualSpaceGetDimension(Qref, &fpdim);
1232: for (i = 0; i < fpdim; ++i) {
1233: PetscDualSpaceGetFunctional(Qref, i, &f);
1234: PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1235: npoints += Np;
1236: }
1237: PetscMalloc1(npoints*dim,&points);
1238: for (i = 0, k = 0; i < fpdim; ++i) {
1239: PetscDualSpaceGetFunctional(Qref, i, &f);
1240: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1241: for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1242: }
1244: for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1245: PetscObject obj;
1246: PetscClassId id;
1247: PetscReal *B;
1248: PetscInt NcJ = 0, cpdim = 0, j, qNc;
1250: PetscDSGetDiscretization(prob, fieldJ, &obj);
1251: PetscObjectGetClassId(obj, &id);
1252: if (id == PETSCFE_CLASSID) {
1253: PetscFE fe = (PetscFE) obj;
1255: /* Evaluate basis at points */
1256: PetscFEGetNumComponents(fe, &NcJ);
1257: PetscFEGetDimension(fe, &cpdim);
1258: /* For now, fields only interpolate themselves */
1259: if (fieldI == fieldJ) {
1260: 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);
1261: PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1262: for (i = 0, k = 0; i < fpdim; ++i) {
1263: PetscDualSpaceGetFunctional(Qref, i, &f);
1264: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1265: 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);
1266: for (p = 0; p < Np; ++p, ++k) {
1267: for (j = 0; j < cpdim; ++j) {
1268: /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */
1269: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1270: }
1271: }
1272: }
1273: PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1274: }
1275: } else if (id == PETSCFV_CLASSID) {
1276: PetscFV fv = (PetscFV) obj;
1278: /* Evaluate constant function at points */
1279: PetscFVGetNumComponents(fv, &NcJ);
1280: cpdim = 1;
1281: /* For now, fields only interpolate themselves */
1282: if (fieldI == fieldJ) {
1283: 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);
1284: for (i = 0, k = 0; i < fpdim; ++i) {
1285: PetscDualSpaceGetFunctional(Qref, i, &f);
1286: PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1287: 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);
1288: for (p = 0; p < Np; ++p, ++k) {
1289: for (j = 0; j < cpdim; ++j) {
1290: for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c];
1291: }
1292: }
1293: }
1294: }
1295: }
1296: offsetJ += cpdim*NcJ;
1297: }
1298: offsetI += fpdim*Nc;
1299: PetscFree(points);
1300: }
1301: if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1302: /* Preallocate matrix */
1303: {
1304: Mat preallocator;
1305: PetscScalar *vals;
1306: PetscInt *cellCIndices, *cellFIndices;
1307: PetscInt locRows, locCols, cell;
1309: MatGetLocalSize(In, &locRows, &locCols);
1310: MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1311: MatSetType(preallocator, MATPREALLOCATOR);
1312: MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1313: MatSetUp(preallocator);
1314: PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1315: for (cell = cStart; cell < cEnd; ++cell) {
1316: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1317: MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1318: }
1319: PetscFree3(vals,cellCIndices,cellFIndices);
1320: MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1321: MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1322: MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1323: MatDestroy(&preallocator);
1324: }
1325: /* Fill matrix */
1326: MatZeroEntries(In);
1327: for (c = cStart; c < cEnd; ++c) {
1328: DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1329: }
1330: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1331: PetscFree2(feRef,fvRef);
1332: PetscFree(elemMat);
1333: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1334: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1335: if (mesh->printFEM) {
1336: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1337: MatChop(In, 1.0e-10);
1338: MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1339: }
1340: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1341: return(0);
1342: }
1344: /*@
1345: DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.
1347: Input Parameters:
1348: + dmf - The fine mesh
1349: . dmc - The coarse mesh
1350: - user - The user context
1352: Output Parameter:
1353: . In - The interpolation matrix
1355: Level: developer
1357: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1358: @*/
1359: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1360: {
1361: DM_Plex *mesh = (DM_Plex *) dmf->data;
1362: const char *name = "Interpolator";
1363: PetscDS prob;
1364: PetscSection fsection, csection, globalFSection, globalCSection;
1365: PetscHashJK ht;
1366: PetscLayout rLayout;
1367: PetscInt *dnz, *onz;
1368: PetscInt locRows, rStart, rEnd;
1369: PetscReal *x, *v0, *J, *invJ, detJ;
1370: PetscReal *v0c, *Jc, *invJc, detJc;
1371: PetscScalar *elemMat;
1372: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell;
1376: PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1377: DMGetCoordinateDim(dmc, &dim);
1378: DMGetDS(dmc, &prob);
1379: PetscDSGetRefCoordArrays(prob, &x, NULL);
1380: PetscDSGetNumFields(prob, &Nf);
1381: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1382: PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1383: DMGetDefaultSection(dmf, &fsection);
1384: DMGetDefaultGlobalSection(dmf, &globalFSection);
1385: DMGetDefaultSection(dmc, &csection);
1386: DMGetDefaultGlobalSection(dmc, &globalCSection);
1387: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1388: PetscDSGetTotalDimension(prob, &totDim);
1389: PetscMalloc1(totDim, &elemMat);
1391: MatGetLocalSize(In, &locRows, NULL);
1392: PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1393: PetscLayoutSetLocalSize(rLayout, locRows);
1394: PetscLayoutSetBlockSize(rLayout, 1);
1395: PetscLayoutSetUp(rLayout);
1396: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1397: PetscLayoutDestroy(&rLayout);
1398: PetscCalloc2(locRows,&dnz,locRows,&onz);
1399: PetscHashJKCreate(&ht);
1400: for (field = 0; field < Nf; ++field) {
1401: PetscObject obj;
1402: PetscClassId id;
1403: PetscDualSpace Q = NULL;
1404: PetscQuadrature f;
1405: const PetscReal *qpoints;
1406: PetscInt Nc, Np, fpdim, i, d;
1408: PetscDSGetDiscretization(prob, field, &obj);
1409: PetscObjectGetClassId(obj, &id);
1410: if (id == PETSCFE_CLASSID) {
1411: PetscFE fe = (PetscFE) obj;
1413: PetscFEGetDualSpace(fe, &Q);
1414: PetscFEGetNumComponents(fe, &Nc);
1415: } else if (id == PETSCFV_CLASSID) {
1416: PetscFV fv = (PetscFV) obj;
1418: PetscFVGetDualSpace(fv, &Q);
1419: Nc = 1;
1420: }
1421: PetscDualSpaceGetDimension(Q, &fpdim);
1422: /* For each fine grid cell */
1423: for (cell = cStart; cell < cEnd; ++cell) {
1424: PetscInt *findices, *cindices;
1425: PetscInt numFIndices, numCIndices;
1427: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1428: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1429: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1430: for (i = 0; i < fpdim; ++i) {
1431: Vec pointVec;
1432: PetscScalar *pV;
1433: PetscSF coarseCellSF = NULL;
1434: const PetscSFNode *coarseCells;
1435: PetscInt numCoarseCells, q, c;
1437: /* Get points from the dual basis functional quadrature */
1438: PetscDualSpaceGetFunctional(Q, i, &f);
1439: PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1440: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1441: VecSetBlockSize(pointVec, dim);
1442: VecGetArray(pointVec, &pV);
1443: for (q = 0; q < Np; ++q) {
1444: /* Transform point to real space */
1445: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1446: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1447: }
1448: VecRestoreArray(pointVec, &pV);
1449: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1450: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1451: PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1452: /* Update preallocation info */
1453: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1454: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1455: {
1456: PetscHashJKKey key;
1457: PetscHashJKIter missing, iter;
1459: key.j = findices[i];
1460: if (key.j >= 0) {
1461: /* Get indices for coarse elements */
1462: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1463: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1464: for (c = 0; c < numCIndices; ++c) {
1465: key.k = cindices[c];
1466: if (key.k < 0) continue;
1467: PetscHashJKPut(ht, key, &missing, &iter);
1468: if (missing) {
1469: PetscHashJKSet(ht, iter, 1);
1470: if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1471: else ++onz[key.j-rStart];
1472: }
1473: }
1474: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1475: }
1476: }
1477: }
1478: PetscSFDestroy(&coarseCellSF);
1479: VecDestroy(&pointVec);
1480: }
1481: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1482: }
1483: }
1484: PetscHashJKDestroy(&ht);
1485: MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1486: MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1487: PetscFree2(dnz,onz);
1488: for (field = 0; field < Nf; ++field) {
1489: PetscObject obj;
1490: PetscClassId id;
1491: PetscDualSpace Q = NULL;
1492: PetscQuadrature f;
1493: const PetscReal *qpoints, *qweights;
1494: PetscInt Nc, qNc, Np, fpdim, i, d;
1496: PetscDSGetDiscretization(prob, field, &obj);
1497: PetscObjectGetClassId(obj, &id);
1498: if (id == PETSCFE_CLASSID) {
1499: PetscFE fe = (PetscFE) obj;
1501: PetscFEGetDualSpace(fe, &Q);
1502: PetscFEGetNumComponents(fe, &Nc);
1503: } else if (id == PETSCFV_CLASSID) {
1504: PetscFV fv = (PetscFV) obj;
1506: PetscFVGetDualSpace(fv, &Q);
1507: Nc = 1;
1508: }
1509: PetscDualSpaceGetDimension(Q, &fpdim);
1510: /* For each fine grid cell */
1511: for (cell = cStart; cell < cEnd; ++cell) {
1512: PetscInt *findices, *cindices;
1513: PetscInt numFIndices, numCIndices;
1515: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1516: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1517: if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1518: for (i = 0; i < fpdim; ++i) {
1519: Vec pointVec;
1520: PetscScalar *pV;
1521: PetscSF coarseCellSF = NULL;
1522: const PetscSFNode *coarseCells;
1523: PetscInt numCoarseCells, cpdim, q, c, j;
1525: /* Get points from the dual basis functional quadrature */
1526: PetscDualSpaceGetFunctional(Q, i, &f);
1527: PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
1528: 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);
1529: VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1530: VecSetBlockSize(pointVec, dim);
1531: VecGetArray(pointVec, &pV);
1532: for (q = 0; q < Np; ++q) {
1533: /* Transform point to real space */
1534: CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1535: for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1536: }
1537: VecRestoreArray(pointVec, &pV);
1538: /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1539: DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1540: /* Update preallocation info */
1541: PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1542: if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1543: VecGetArray(pointVec, &pV);
1544: for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1545: PetscReal pVReal[3];
1547: DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1548: /* Transform points from real space to coarse reference space */
1549: DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
1550: for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
1551: CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x);
1553: if (id == PETSCFE_CLASSID) {
1554: PetscFE fe = (PetscFE) obj;
1555: PetscReal *B;
1557: /* Evaluate coarse basis on contained point */
1558: PetscFEGetDimension(fe, &cpdim);
1559: PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1560: PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
1561: /* Get elemMat entries by multiplying by weight */
1562: for (j = 0; j < cpdim; ++j) {
1563: for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
1564: }
1565: PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1566: } else {
1567: cpdim = 1;
1568: for (j = 0; j < cpdim; ++j) {
1569: for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
1570: }
1571: }
1572: /* Update interpolator */
1573: if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
1574: if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
1575: MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
1576: DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1577: }
1578: VecRestoreArray(pointVec, &pV);
1579: PetscSFDestroy(&coarseCellSF);
1580: VecDestroy(&pointVec);
1581: }
1582: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1583: }
1584: }
1585: PetscFree3(v0,J,invJ);
1586: PetscFree3(v0c,Jc,invJc);
1587: PetscFree(elemMat);
1588: MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1589: MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1590: PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1591: return(0);
1592: }
1594: /*@
1595: DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns
1597: Input Parameters:
1598: + dmc - The coarse mesh
1599: - dmf - The fine mesh
1600: - user - The user context
1602: Output Parameter:
1603: . sc - The mapping
1605: Level: developer
1607: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1608: @*/
1609: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1610: {
1611: PetscDS prob;
1612: PetscFE *feRef;
1613: PetscFV *fvRef;
1614: Vec fv, cv;
1615: IS fis, cis;
1616: PetscSection fsection, fglobalSection, csection, cglobalSection;
1617: PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1618: PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
1622: PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1623: DMGetDimension(dmf, &dim);
1624: DMGetDefaultSection(dmf, &fsection);
1625: DMGetDefaultGlobalSection(dmf, &fglobalSection);
1626: DMGetDefaultSection(dmc, &csection);
1627: DMGetDefaultGlobalSection(dmc, &cglobalSection);
1628: PetscSectionGetNumFields(fsection, &Nf);
1629: DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1630: DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1631: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1632: DMGetDS(dmc, &prob);
1633: PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1634: for (f = 0; f < Nf; ++f) {
1635: PetscObject obj;
1636: PetscClassId id;
1637: PetscInt fNb = 0, Nc = 0;
1639: PetscDSGetDiscretization(prob, f, &obj);
1640: PetscObjectGetClassId(obj, &id);
1641: if (id == PETSCFE_CLASSID) {
1642: PetscFE fe = (PetscFE) obj;
1644: PetscFERefine(fe, &feRef[f]);
1645: PetscFEGetDimension(feRef[f], &fNb);
1646: PetscFEGetNumComponents(fe, &Nc);
1647: } else if (id == PETSCFV_CLASSID) {
1648: PetscFV fv = (PetscFV) obj;
1649: PetscDualSpace Q;
1651: PetscFVRefine(fv, &fvRef[f]);
1652: PetscFVGetDualSpace(fvRef[f], &Q);
1653: PetscDualSpaceGetDimension(Q, &fNb);
1654: PetscFVGetNumComponents(fv, &Nc);
1655: }
1656: fTotDim += fNb*Nc;
1657: }
1658: PetscDSGetTotalDimension(prob, &cTotDim);
1659: PetscMalloc1(cTotDim,&cmap);
1660: for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1661: PetscFE feC;
1662: PetscFV fvC;
1663: PetscDualSpace QF, QC;
1664: PetscInt NcF, NcC, fpdim, cpdim;
1666: if (feRef[field]) {
1667: PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1668: PetscFEGetNumComponents(feC, &NcC);
1669: PetscFEGetNumComponents(feRef[field], &NcF);
1670: PetscFEGetDualSpace(feRef[field], &QF);
1671: PetscDualSpaceGetDimension(QF, &fpdim);
1672: PetscFEGetDualSpace(feC, &QC);
1673: PetscDualSpaceGetDimension(QC, &cpdim);
1674: } else {
1675: PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1676: PetscFVGetNumComponents(fvC, &NcC);
1677: PetscFVGetNumComponents(fvRef[field], &NcF);
1678: PetscFVGetDualSpace(fvRef[field], &QF);
1679: PetscDualSpaceGetDimension(QF, &fpdim);
1680: PetscFVGetDualSpace(fvC, &QC);
1681: PetscDualSpaceGetDimension(QC, &cpdim);
1682: }
1683: 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);
1684: for (c = 0; c < cpdim; ++c) {
1685: PetscQuadrature cfunc;
1686: const PetscReal *cqpoints;
1687: PetscInt NpC;
1688: PetscBool found = PETSC_FALSE;
1690: PetscDualSpaceGetFunctional(QC, c, &cfunc);
1691: PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);
1692: if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1693: for (f = 0; f < fpdim; ++f) {
1694: PetscQuadrature ffunc;
1695: const PetscReal *fqpoints;
1696: PetscReal sum = 0.0;
1697: PetscInt NpF, comp;
1699: PetscDualSpaceGetFunctional(QF, f, &ffunc);
1700: PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);
1701: if (NpC != NpF) continue;
1702: for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1703: if (sum > 1.0e-9) continue;
1704: for (comp = 0; comp < NcC; ++comp) {
1705: cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1706: }
1707: found = PETSC_TRUE;
1708: break;
1709: }
1710: if (!found) {
1711: /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1712: if (fvRef[field]) {
1713: PetscInt comp;
1714: for (comp = 0; comp < NcC; ++comp) {
1715: cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1716: }
1717: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1718: }
1719: }
1720: offsetC += cpdim*NcC;
1721: offsetF += fpdim*NcF;
1722: }
1723: for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1724: PetscFree2(feRef,fvRef);
1726: DMGetGlobalVector(dmf, &fv);
1727: DMGetGlobalVector(dmc, &cv);
1728: VecGetOwnershipRange(cv, &startC, &endC);
1729: PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1730: PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1731: PetscMalloc1(m,&cindices);
1732: PetscMalloc1(m,&findices);
1733: for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1734: for (c = cStart; c < cEnd; ++c) {
1735: DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1736: for (d = 0; d < cTotDim; ++d) {
1737: if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
1738: 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]]);
1739: cindices[cellCIndices[d]-startC] = cellCIndices[d];
1740: findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1741: }
1742: }
1743: PetscFree(cmap);
1744: PetscFree2(cellCIndices,cellFIndices);
1746: ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1747: ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1748: VecScatterCreate(cv, cis, fv, fis, sc);
1749: ISDestroy(&cis);
1750: ISDestroy(&fis);
1751: DMRestoreGlobalVector(dmf, &fv);
1752: DMRestoreGlobalVector(dmc, &cv);
1753: PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1754: return(0);
1755: }