Actual source code: dmceed.c
1: #include <petsc/private/dmimpl.h>
2: #include <petscdmceed.h>
4: #ifdef PETSC_HAVE_LIBCEED
5: #include <petsc/private/dmpleximpl.h>
6: #include <petscdmplexceed.h>
7: #include <petscfeceed.h>
9: /*@C
10: DMGetCeed - Get the LibCEED context associated with this `DM`
12: Not Collective
14: Input Parameter:
15: . DM - The `DM`
17: Output Parameter:
18: . ceed - The LibCEED context
20: Level: intermediate
22: .seealso: `DM`, `DMCreate()`
23: @*/
24: PetscErrorCode DMGetCeed(DM dm, Ceed *ceed)
25: {
26: PetscFunctionBegin;
28: PetscAssertPointer(ceed, 2);
29: if (!dm->ceed) {
30: char ceedresource[PETSC_MAX_PATH_LEN]; /* libCEED resource specifier */
31: const char *prefix;
33: PetscCall(PetscStrncpy(ceedresource, "/cpu/self", sizeof(ceedresource)));
34: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
35: PetscCall(PetscOptionsGetString(NULL, prefix, "-dm_ceed", ceedresource, sizeof(ceedresource), NULL));
36: PetscCallCEED(CeedInit(ceedresource, &dm->ceed));
37: }
38: *ceed = dm->ceed;
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static CeedMemType PetscMemType2Ceed(PetscMemType mem_type)
43: {
44: return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST;
45: }
47: PetscErrorCode VecGetCeedVector(Vec X, Ceed ceed, CeedVector *cx)
48: {
49: PetscMemType memtype;
50: PetscScalar *x;
51: PetscInt n;
53: PetscFunctionBegin;
54: PetscCall(VecGetLocalSize(X, &n));
55: PetscCall(VecGetArrayAndMemType(X, &x, &memtype));
56: PetscCallCEED(CeedVectorCreate(ceed, n, cx));
57: PetscCallCEED(CeedVectorSetArray(*cx, PetscMemType2Ceed(memtype), CEED_USE_POINTER, x));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: PetscErrorCode VecRestoreCeedVector(Vec X, CeedVector *cx)
62: {
63: PetscFunctionBegin;
64: PetscCall(VecRestoreArrayAndMemType(X, NULL));
65: PetscCallCEED(CeedVectorDestroy(cx));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: PetscErrorCode VecGetCeedVectorRead(Vec X, Ceed ceed, CeedVector *cx)
70: {
71: PetscMemType memtype;
72: const PetscScalar *x;
73: PetscInt n;
75: PetscFunctionBegin;
76: PetscCall(VecGetLocalSize(X, &n));
77: PetscCall(VecGetArrayReadAndMemType(X, &x, &memtype));
78: PetscCallCEED(CeedVectorCreate(ceed, n, cx));
79: PetscCallCEED(CeedVectorSetArray(*cx, PetscMemType2Ceed(memtype), CEED_USE_POINTER, (PetscScalar *)x));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: PetscErrorCode VecRestoreCeedVectorRead(Vec X, CeedVector *cx)
84: {
85: PetscFunctionBegin;
86: PetscCall(VecRestoreArrayReadAndMemType(X, NULL));
87: PetscCallCEED(CeedVectorDestroy(cx));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: CEED_QFUNCTION(Geometry2D)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
92: {
93: const CeedScalar *x = in[0], *Jac = in[1], *w = in[2];
94: CeedScalar *qdata = out[0];
96: CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
97: {
98: const CeedScalar J[2][2] = {
99: {Jac[i + Q * 0], Jac[i + Q * 2]},
100: {Jac[i + Q * 1], Jac[i + Q * 3]}
101: };
102: const CeedScalar det = J[0][0] * J[1][1] - J[0][1] * J[1][0];
104: qdata[i + Q * 0] = det * w[i];
105: qdata[i + Q * 1] = x[i + Q * 0];
106: qdata[i + Q * 2] = x[i + Q * 1];
107: qdata[i + Q * 3] = J[1][1] / det;
108: qdata[i + Q * 4] = -J[1][0] / det;
109: qdata[i + Q * 5] = -J[0][1] / det;
110: qdata[i + Q * 6] = J[0][0] / det;
111: }
112: return CEED_ERROR_SUCCESS;
113: }
115: CEED_QFUNCTION(Geometry3D)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
116: {
117: const CeedScalar *Jac = in[1], *w = in[2];
118: CeedScalar *qdata = out[0];
120: CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
121: {
122: const CeedScalar J[3][3] = {
123: {Jac[i + Q * 0], Jac[i + Q * 3], Jac[i + Q * 6]},
124: {Jac[i + Q * 1], Jac[i + Q * 4], Jac[i + Q * 7]},
125: {Jac[i + Q * 2], Jac[i + Q * 5], Jac[i + Q * 8]}
126: };
127: const CeedScalar det = J[0][0] * (J[1][1] * J[2][2] - J[1][2] * J[2][1]) + J[0][1] * (J[1][2] * J[2][0] - J[1][0] * J[2][2]) + J[0][2] * (J[1][0] * J[2][1] - J[1][1] * J[2][0]);
129: qdata[i + Q * 0] = det * w[i]; /* det J * weight */
130: }
131: return CEED_ERROR_SUCCESS;
132: }
134: static PetscErrorCode DMCeedCreateGeometry(DM dm, IS cellIS, PetscInt *Nqdata, CeedElemRestriction *erq, CeedVector *qd, DMCeed *soldata)
135: {
136: Ceed ceed;
137: DMCeed sd;
138: PetscDS ds;
139: PetscFE fe;
140: CeedQFunctionUser geom = NULL;
141: const char *geomName = NULL;
142: const PetscInt *cells;
143: PetscInt dim, cdim, cStart, cEnd, Ncell;
144: CeedInt Nq;
146: PetscFunctionBegin;
147: PetscCall(PetscCalloc1(1, &sd));
148: PetscCall(DMGetDimension(dm, &dim));
149: PetscCall(DMGetCoordinateDim(dm, &cdim));
150: PetscCall(DMGetCeed(dm, &ceed));
151: PetscCall(ISGetPointRange(cellIS, &cStart, &cEnd, &cells));
152: Ncell = cEnd - cStart;
154: PetscCall(DMGetCellDS(dm, cells ? cells[cStart] : cStart, &ds, NULL));
155: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
156: PetscCall(PetscFEGetCeedBasis(fe, &sd->basis));
157: PetscCall(CeedBasisGetNumQuadraturePoints(sd->basis, &Nq));
158: PetscCall(DMPlexGetCeedRestriction(dm, NULL, 0, 0, 0, &sd->er));
160: *Nqdata = 1 + cdim + cdim * dim; // |J| * w_q, x, J^{-1}
161: PetscCallCEED(CeedElemRestrictionCreateStrided(ceed, Ncell, Nq, *Nqdata, Ncell * Nq * (*Nqdata), CEED_STRIDES_BACKEND, erq));
163: switch (dim) {
164: case 2:
165: geom = Geometry2D;
166: geomName = Geometry2D_loc;
167: break;
168: case 3:
169: geom = Geometry3D;
170: geomName = Geometry3D_loc;
171: break;
172: }
173: PetscCallCEED(CeedQFunctionCreateInterior(ceed, 1, geom, geomName, &sd->qf));
174: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "x", cdim, CEED_EVAL_INTERP));
175: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "dx", cdim * dim, CEED_EVAL_GRAD));
176: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "weight", 1, CEED_EVAL_WEIGHT));
177: PetscCallCEED(CeedQFunctionAddOutput(sd->qf, "qdata", *Nqdata, CEED_EVAL_NONE));
179: PetscCallCEED(CeedOperatorCreate(ceed, sd->qf, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &sd->op));
180: PetscCallCEED(CeedOperatorSetField(sd->op, "x", sd->er, sd->basis, CEED_VECTOR_ACTIVE));
181: PetscCallCEED(CeedOperatorSetField(sd->op, "dx", sd->er, sd->basis, CEED_VECTOR_ACTIVE));
182: PetscCallCEED(CeedOperatorSetField(sd->op, "weight", CEED_ELEMRESTRICTION_NONE, sd->basis, CEED_VECTOR_NONE));
183: PetscCallCEED(CeedOperatorSetField(sd->op, "qdata", *erq, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
185: PetscCallCEED(CeedElemRestrictionCreateVector(*erq, qd, NULL));
186: *soldata = sd;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: PetscErrorCode DMRefineHook_Ceed(DM coarse, DM fine, void *ctx)
191: {
192: PetscFunctionBegin;
193: if (coarse->dmceed) PetscCall(DMCeedCreate(fine, coarse->dmceed->geom ? PETSC_TRUE : PETSC_FALSE, coarse->dmceed->func, coarse->dmceed->funcSource));
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: PetscErrorCode DMCeedCreate_Internal(DM dm, IS cellIS, PetscBool createGeometry, CeedQFunctionUser func, const char *func_source, DMCeed *soldata)
198: {
199: PetscDS ds;
200: PetscFE fe;
201: DMCeed sd;
202: Ceed ceed;
203: PetscInt dim, Nc, Nqdata = 0;
204: CeedInt Nq;
206: PetscFunctionBegin;
207: PetscCall(PetscCalloc1(1, &sd));
208: PetscCall(DMGetDimension(dm, &dim));
209: PetscCall(DMGetCeed(dm, &ceed));
210: PetscCall(DMGetDS(dm, &ds));
211: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe));
212: PetscCall(PetscFEGetCeedBasis(fe, &sd->basis));
213: PetscCall(PetscFEGetNumComponents(fe, &Nc));
214: PetscCall(CeedBasisGetNumQuadraturePoints(sd->basis, &Nq));
215: PetscCall(DMPlexGetCeedRestriction(dm, NULL, 0, 0, 0, &sd->er));
217: if (createGeometry) {
218: DM cdm;
220: PetscCall(DMGetCoordinateDM(dm, &cdm));
221: PetscCall(DMCeedCreateGeometry(cdm, cellIS, &Nqdata, &sd->erq, &sd->qd, &sd->geom));
222: }
224: if (sd->geom) {
225: PetscInt cdim;
226: CeedInt Nqx;
228: PetscCallCEED(CeedBasisGetNumQuadraturePoints(sd->geom->basis, &Nqx));
229: PetscCheck(Nqx == Nq, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP, "Number of qpoints for solution %" CeedInt_FMT " != %" CeedInt_FMT " Number of qpoints for coordinates", Nq, Nqx);
230: /* TODO Remove this limitation */
231: PetscCall(DMGetCoordinateDim(dm, &cdim));
232: PetscCheck(dim == cdim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_INCOMP, "Topological dimension %" PetscInt_FMT " != %" PetscInt_FMT " embedding dimension", dim, cdim);
233: }
235: PetscCallCEED(CeedQFunctionCreateInterior(ceed, 1, func, func_source, &sd->qf));
236: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "u", Nc, CEED_EVAL_INTERP));
237: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "du", Nc * dim, CEED_EVAL_GRAD));
238: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "qdata", Nqdata, CEED_EVAL_NONE));
239: PetscCallCEED(CeedQFunctionAddOutput(sd->qf, "v", Nc, CEED_EVAL_INTERP));
240: PetscCallCEED(CeedQFunctionAddOutput(sd->qf, "dv", Nc * dim, CEED_EVAL_GRAD));
242: PetscCallCEED(CeedOperatorCreate(ceed, sd->qf, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &sd->op));
243: PetscCallCEED(CeedOperatorSetField(sd->op, "u", sd->er, sd->basis, CEED_VECTOR_ACTIVE));
244: PetscCallCEED(CeedOperatorSetField(sd->op, "du", sd->er, sd->basis, CEED_VECTOR_ACTIVE));
245: PetscCallCEED(CeedOperatorSetField(sd->op, "qdata", sd->erq, CEED_BASIS_NONE, sd->qd));
246: PetscCallCEED(CeedOperatorSetField(sd->op, "v", sd->er, sd->basis, CEED_VECTOR_ACTIVE));
247: PetscCallCEED(CeedOperatorSetField(sd->op, "dv", sd->er, sd->basis, CEED_VECTOR_ACTIVE));
249: // Handle refinement
250: sd->func = func;
251: PetscCall(PetscStrallocpy(func_source, &sd->funcSource));
252: PetscCall(DMRefineHookAdd(dm, DMRefineHook_Ceed, NULL, NULL));
254: *soldata = sd;
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: PetscErrorCode DMCeedCreate(DM dm, PetscBool createGeometry, CeedQFunctionUser func, const char *func_source)
259: {
260: DM plex;
261: IS cellIS;
263: PetscFunctionBegin;
264: PetscCall(DMConvert(dm, DMPLEX, &plex));
265: PetscCall(DMPlexGetAllCells_Internal(plex, &cellIS));
266: #ifdef PETSC_HAVE_LIBCEED
267: PetscCall(DMCeedCreate_Internal(dm, cellIS, createGeometry, func, func_source, &dm->dmceed));
268: #endif
269: PetscCall(ISDestroy(&cellIS));
270: PetscCall(DMDestroy(&plex));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: static PetscErrorCode DMCeedCreateGeometryFVM(DM dm, IS faceIS, PetscInt *Nqdata, CeedElemRestriction *erq, CeedVector *qd, DMCeed *soldata)
275: {
276: Ceed ceed;
277: DMCeed sd;
278: const PetscInt *faces;
279: CeedInt strides[3];
280: PetscInt dim, cdim, fStart, fEnd, Nface, Nq = 1;
282: PetscFunctionBegin;
283: PetscCall(PetscCalloc1(1, &sd));
284: PetscCall(DMGetDimension(dm, &dim));
285: PetscCall(DMGetCoordinateDim(dm, &cdim));
286: PetscCall(DMGetCeed(dm, &ceed));
287: PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
288: Nface = fEnd - fStart;
290: *Nqdata = cdim + 2; // face normal and support cell volumes
291: strides[0] = 1;
292: strides[1] = Nq;
293: strides[2] = Nq * (*Nqdata);
294: PetscCallCEED(CeedElemRestrictionCreateStrided(ceed, Nface, Nq, *Nqdata, Nface * Nq * (*Nqdata), strides, erq));
295: PetscCallCEED(CeedElemRestrictionCreateVector(*erq, qd, NULL));
296: *soldata = sd;
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode DMCeedCreateInfoFVM(DM dm, IS faceIS, PetscInt *Nqinfo, CeedElemRestriction *eri, CeedVector *qi, DMCeed *solinfo)
301: {
302: Ceed ceed;
303: DMCeed si;
304: const PetscInt *faces;
305: CeedInt strides[3];
306: PetscInt dim, fStart, fEnd, Nface, Nq = 1;
308: PetscFunctionBegin;
309: PetscCall(PetscCalloc1(1, &si));
310: PetscCall(DMGetDimension(dm, &dim));
311: PetscCall(DMGetCeed(dm, &ceed));
312: PetscCall(ISGetPointRange(faceIS, &fStart, &fEnd, &faces));
313: Nface = fEnd - fStart;
315: *Nqinfo = 3; // face number and support cell numbers
316: strides[0] = 1;
317: strides[1] = Nq;
318: strides[2] = Nq * (*Nqinfo);
319: PetscCallCEED(CeedElemRestrictionCreateStrided(ceed, Nface, Nq, *Nqinfo, Nface * Nq * (*Nqinfo), strides, eri));
320: PetscCallCEED(CeedElemRestrictionCreateVector(*eri, qi, NULL));
321: *solinfo = si;
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: PetscErrorCode DMCeedCreateFVM_Internal(DM dm, IS faceIS, PetscBool createGeometry, PetscBool createInfo, CeedQFunctionUser func, const char *func_source, DMCeed *soldata, CeedQFunctionContext qfCtx)
326: {
327: PetscDS ds;
328: PetscFV fv;
329: DMCeed sd;
330: Ceed ceed;
331: PetscInt dim, Nc, Nqdata = 0, Nqinfo = 0;
333: PetscFunctionBegin;
334: PetscCall(PetscCalloc1(1, &sd));
335: PetscCall(DMGetDimension(dm, &dim));
336: PetscCall(DMGetCeed(dm, &ceed));
337: PetscCall(DMGetDS(dm, &ds));
338: PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fv));
339: PetscCall(PetscFVGetNumComponents(fv, &Nc));
340: PetscCall(DMPlexCreateCeedRestrictionFVM(dm, &sd->erL, &sd->erR));
342: if (createGeometry) {
343: DM cdm;
345: PetscCall(DMGetCoordinateDM(dm, &cdm));
346: PetscCall(DMCeedCreateGeometryFVM(cdm, faceIS, &Nqdata, &sd->erq, &sd->qd, &sd->geom));
347: }
349: if (createInfo) {
350: DM cdm;
352: PetscCall(DMGetCoordinateDM(dm, &cdm));
353: PetscCall(DMCeedCreateInfoFVM(cdm, faceIS, &Nqinfo, &sd->eri, &sd->qi, &sd->info));
354: PetscCall(DMCeedComputeInfo(dm, sd));
355: }
357: PetscCallCEED(CeedQFunctionCreateInterior(ceed, 1, func, func_source, &sd->qf));
358: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "uL", Nc, CEED_EVAL_NONE));
359: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "uR", Nc, CEED_EVAL_NONE));
360: PetscCallCEED(CeedQFunctionAddInput(sd->qf, "geom", Nqdata, CEED_EVAL_NONE));
361: if (createInfo) PetscCallCEED(CeedQFunctionAddInput(sd->qf, "info", Nqinfo, CEED_EVAL_NONE));
362: PetscCallCEED(CeedQFunctionAddOutput(sd->qf, "cL", Nc, CEED_EVAL_NONE));
363: PetscCallCEED(CeedQFunctionAddOutput(sd->qf, "cR", Nc, CEED_EVAL_NONE));
365: PetscCallCEED(CeedQFunctionSetContext(sd->qf, qfCtx));
367: PetscCallCEED(CeedOperatorCreate(ceed, sd->qf, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &sd->op));
368: PetscCallCEED(CeedOperatorSetField(sd->op, "uL", sd->erL, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
369: PetscCallCEED(CeedOperatorSetField(sd->op, "uR", sd->erR, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
370: PetscCallCEED(CeedOperatorSetField(sd->op, "geom", sd->erq, CEED_BASIS_NONE, sd->qd));
371: if (createInfo) PetscCallCEED(CeedOperatorSetField(sd->op, "info", sd->eri, CEED_BASIS_NONE, sd->qi));
372: PetscCallCEED(CeedOperatorSetField(sd->op, "cL", sd->erL, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
373: PetscCallCEED(CeedOperatorSetField(sd->op, "cR", sd->erR, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
375: // Handle refinement
376: sd->func = func;
377: PetscCall(PetscStrallocpy(func_source, &sd->funcSource));
378: PetscCall(DMRefineHookAdd(dm, DMRefineHook_Ceed, NULL, NULL));
380: *soldata = sd;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: PetscErrorCode DMCeedCreateFVM(DM dm, PetscBool createGeometry, CeedQFunctionUser func, const char *func_source, CeedQFunctionContext qfCtx)
385: {
386: DM plex;
387: IS faceIS;
389: PetscFunctionBegin;
390: PetscCall(DMConvert(dm, DMPLEX, &plex));
391: PetscCall(DMPlexGetAllFaces_Internal(plex, &faceIS));
392: #ifdef PETSC_HAVE_LIBCEED
393: PetscCall(DMCeedCreateFVM_Internal(dm, faceIS, createGeometry, PETSC_TRUE, func, func_source, &dm->dmceed, qfCtx));
394: #endif
395: PetscCall(ISDestroy(&faceIS));
396: PetscCall(DMDestroy(&plex));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: #endif
402: PetscErrorCode DMCeedDestroy(DMCeed *pceed)
403: {
404: DMCeed p = *pceed;
406: PetscFunctionBegin;
407: if (!p) PetscFunctionReturn(PETSC_SUCCESS);
408: #ifdef PETSC_HAVE_LIBCEED
409: PetscCall(PetscFree(p->funcSource));
410: if (p->qd) PetscCallCEED(CeedVectorDestroy(&p->qd));
411: if (p->qi) PetscCallCEED(CeedVectorDestroy(&p->qi));
412: if (p->op) PetscCallCEED(CeedOperatorDestroy(&p->op));
413: if (p->qf) PetscCallCEED(CeedQFunctionDestroy(&p->qf));
414: if (p->erL) PetscCallCEED(CeedElemRestrictionDestroy(&p->erL));
415: if (p->erR) PetscCallCEED(CeedElemRestrictionDestroy(&p->erR));
416: if (p->erq) PetscCallCEED(CeedElemRestrictionDestroy(&p->erq));
417: if (p->eri) PetscCallCEED(CeedElemRestrictionDestroy(&p->eri));
418: PetscCall(DMCeedDestroy(&p->geom));
419: PetscCall(DMCeedDestroy(&p->info));
420: #endif
421: PetscCall(PetscFree(p));
422: *pceed = NULL;
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: PetscErrorCode DMCeedComputeGeometry(DM dm, DMCeed sd)
427: {
428: #ifdef PETSC_HAVE_LIBCEED
429: Ceed ceed;
430: Vec coords;
431: CeedVector ccoords;
432: #endif
434: PetscFunctionBegin;
435: #ifdef PETSC_HAVE_LIBCEED
436: PetscCall(DMGetCeed(dm, &ceed));
437: PetscCall(DMGetCoordinatesLocal(dm, &coords));
438: PetscCall(VecGetCeedVectorRead(coords, ceed, &ccoords));
439: if (sd->geom->op) PetscCallCEED(CeedOperatorApply(sd->geom->op, ccoords, sd->qd, CEED_REQUEST_IMMEDIATE));
440: else PetscCall(DMPlexCeedComputeGeometryFVM(dm, sd->qd));
441: //PetscCallCEED(CeedVectorView(sd->qd, "%g", stdout));
442: PetscCall(VecRestoreCeedVectorRead(coords, &ccoords));
443: #endif
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: PetscErrorCode DMCeedComputeInfo(DM dm, DMCeed sd)
448: {
449: #ifdef PETSC_HAVE_LIBCEED
450: CeedScalar *a;
451: #endif
453: PetscFunctionBegin;
454: #ifdef PETSC_HAVE_LIBCEED
455: PetscCallCEED(CeedVectorGetArrayWrite(sd->qi, CEED_MEM_HOST, &a));
457: IS iterIS;
458: DMLabel label = NULL;
459: const PetscInt *indices;
460: PetscInt value = 0, height = 1, NfInt = 0, Nf = 0;
462: PetscCall(DMGetPoints_Internal(dm, label, value, height, &iterIS));
463: if (iterIS) {
464: PetscCall(ISGetIndices(iterIS, &indices));
465: PetscCall(ISGetLocalSize(iterIS, &Nf));
466: for (PetscInt p = 0, Ns; p < Nf; ++p) {
467: PetscCall(DMPlexGetSupportSize(dm, indices[p], &Ns));
468: if (Ns == 2) ++NfInt;
469: }
470: } else {
471: indices = NULL;
472: }
474: PetscInt infoOffset = 0;
476: for (PetscInt p = 0; p < Nf; ++p) {
477: const PetscInt face = indices[p];
478: const PetscInt *supp;
479: PetscInt Ns;
481: PetscCall(DMPlexGetSupport(dm, face, &supp));
482: PetscCall(DMPlexGetSupportSize(dm, face, &Ns));
483: // Ignore boundary faces
484: // TODO check for face on parallel boundary
485: if (Ns == 2) {
486: a[infoOffset++] = face;
487: a[infoOffset++] = supp[0];
488: a[infoOffset++] = supp[1];
489: }
490: }
491: PetscCheck(infoOffset == NfInt * 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Shape mismatch, info offsets array of shape (%" PetscInt_FMT ") initialized for %" PetscInt_FMT " nodes", infoOffset, NfInt * 3);
492: if (iterIS) PetscCall(ISRestoreIndices(iterIS, &indices));
493: PetscCall(ISDestroy(&iterIS));
495: PetscCallCEED(CeedVectorRestoreArray(sd->qi, &a));
496: #endif
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }