1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
2 #include <petscdmceed.h>
3
4 #ifdef PETSC_HAVE_LIBCEED
5 #include <petsc/private/dmpleximpl.h>
6 #include <petscdmplexceed.h>
7 #include <petscfeceed.h>
8
9 /*@C
10 DMGetCeed - Get the LibCEED context associated with this `DM`
11
12 Not Collective
13
14 Input Parameter:
15 . DM - The `DM`
16
17 Output Parameter:
18 . ceed - The LibCEED context
19
20 Level: intermediate
21
22 .seealso: `DM`, `DMCreate()`
23 @*/
DMGetCeed(DM dm,Ceed * ceed)24 PetscErrorCode DMGetCeed(DM dm, Ceed *ceed)
25 {
26 PetscFunctionBegin;
27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
28 PetscAssertPointer(ceed, 2);
29 if (!dm->ceed) {
30 char ceedresource[PETSC_MAX_PATH_LEN]; /* libCEED resource specifier */
31 const char *prefix;
32
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 }
41
PetscMemType2Ceed(PetscMemType mem_type)42 static CeedMemType PetscMemType2Ceed(PetscMemType mem_type)
43 {
44 return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST;
45 }
46
VecGetCeedVector(Vec X,Ceed ceed,CeedVector * cx)47 PetscErrorCode VecGetCeedVector(Vec X, Ceed ceed, CeedVector *cx)
48 {
49 PetscMemType memtype;
50 PetscScalar *x;
51 PetscInt n;
52
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 }
60
VecRestoreCeedVector(Vec X,CeedVector * cx)61 PetscErrorCode VecRestoreCeedVector(Vec X, CeedVector *cx)
62 {
63 PetscFunctionBegin;
64 PetscCall(VecRestoreArrayAndMemType(X, NULL));
65 PetscCallCEED(CeedVectorDestroy(cx));
66 PetscFunctionReturn(PETSC_SUCCESS);
67 }
68
VecGetCeedVectorRead(Vec X,Ceed ceed,CeedVector * cx)69 PetscErrorCode VecGetCeedVectorRead(Vec X, Ceed ceed, CeedVector *cx)
70 {
71 PetscMemType memtype;
72 const PetscScalar *x;
73 PetscInt n;
74
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 }
82
VecRestoreCeedVectorRead(Vec X,CeedVector * cx)83 PetscErrorCode VecRestoreCeedVectorRead(Vec X, CeedVector *cx)
84 {
85 PetscFunctionBegin;
86 PetscCall(VecRestoreArrayReadAndMemType(X, NULL));
87 PetscCallCEED(CeedVectorDestroy(cx));
88 PetscFunctionReturn(PETSC_SUCCESS);
89 }
90
Geometry2D(PetscCtx ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)91 CEED_QFUNCTION(Geometry2D)(PetscCtx 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];
95
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];
103
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 }
114
Geometry3D(PetscCtx ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)115 CEED_QFUNCTION(Geometry3D)(PetscCtx 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];
119
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]);
128
129 qdata[i + Q * 0] = det * w[i]; /* det J * weight */
130 }
131 return CEED_ERROR_SUCCESS;
132 }
133
DMCeedCreateGeometry(DM dm,IS cellIS,PetscInt * Nqdata,CeedElemRestriction * erq,CeedVector * qd,DMCeed * soldata)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;
145
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;
153
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));
159
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));
162
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));
178
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));
184
185 PetscCallCEED(CeedElemRestrictionCreateVector(*erq, qd, NULL));
186 *soldata = sd;
187 PetscFunctionReturn(PETSC_SUCCESS);
188 }
189
DMRefineHook_Ceed(DM coarse,DM fine,PetscCtx ctx)190 PetscErrorCode DMRefineHook_Ceed(DM coarse, DM fine, PetscCtx 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 }
196
DMCeedCreate_Internal(DM dm,IS cellIS,PetscBool createGeometry,CeedQFunctionUser func,const char * func_source,DMCeed * soldata)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;
205
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));
216
217 if (createGeometry) {
218 DM cdm;
219
220 PetscCall(DMGetCoordinateDM(dm, &cdm));
221 PetscCall(DMCeedCreateGeometry(cdm, cellIS, &Nqdata, &sd->erq, &sd->qd, &sd->geom));
222 }
223
224 if (sd->geom) {
225 PetscInt cdim;
226 CeedInt Nqx;
227
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 }
234
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));
241
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));
248
249 // Handle refinement
250 sd->func = func;
251 PetscCall(PetscStrallocpy(func_source, &sd->funcSource));
252 PetscCall(DMRefineHookAdd(dm, DMRefineHook_Ceed, NULL, NULL));
253
254 *soldata = sd;
255 PetscFunctionReturn(PETSC_SUCCESS);
256 }
257
DMCeedCreate(DM dm,PetscBool createGeometry,CeedQFunctionUser func,const char * func_source)258 PetscErrorCode DMCeedCreate(DM dm, PetscBool createGeometry, CeedQFunctionUser func, const char *func_source)
259 {
260 DM plex;
261 IS cellIS;
262
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 }
273
DMCeedCreateGeometryFVM(DM dm,IS faceIS,PetscInt * Nqdata,CeedElemRestriction * erq,CeedVector * qd,DMCeed * soldata)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;
281
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;
289
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 }
299
DMCeedCreateInfoFVM(DM dm,IS faceIS,PetscInt * Nqinfo,CeedElemRestriction * eri,CeedVector * qi,DMCeed * solinfo)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;
307
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;
314
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 }
324
DMCeedCreateFVM_Internal(DM dm,IS faceIS,PetscBool createGeometry,PetscBool createInfo,CeedQFunctionUser func,const char * func_source,DMCeed * soldata,CeedQFunctionContext qfCtx)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;
332
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));
341
342 if (createGeometry) {
343 DM cdm;
344
345 PetscCall(DMGetCoordinateDM(dm, &cdm));
346 PetscCall(DMCeedCreateGeometryFVM(cdm, faceIS, &Nqdata, &sd->erq, &sd->qd, &sd->geom));
347 }
348
349 if (createInfo) {
350 DM cdm;
351
352 PetscCall(DMGetCoordinateDM(dm, &cdm));
353 PetscCall(DMCeedCreateInfoFVM(cdm, faceIS, &Nqinfo, &sd->eri, &sd->qi, &sd->info));
354 PetscCall(DMCeedComputeInfo(dm, sd));
355 }
356
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));
364
365 PetscCallCEED(CeedQFunctionSetContext(sd->qf, qfCtx));
366
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));
374
375 // Handle refinement
376 sd->func = func;
377 PetscCall(PetscStrallocpy(func_source, &sd->funcSource));
378 PetscCall(DMRefineHookAdd(dm, DMRefineHook_Ceed, NULL, NULL));
379
380 *soldata = sd;
381 PetscFunctionReturn(PETSC_SUCCESS);
382 }
383
DMCeedCreateFVM(DM dm,PetscBool createGeometry,CeedQFunctionUser func,const char * func_source,CeedQFunctionContext qfCtx)384 PetscErrorCode DMCeedCreateFVM(DM dm, PetscBool createGeometry, CeedQFunctionUser func, const char *func_source, CeedQFunctionContext qfCtx)
385 {
386 DM plex;
387 IS faceIS;
388
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 }
399
400 #endif
401
DMCeedDestroy(DMCeed * pceed)402 PetscErrorCode DMCeedDestroy(DMCeed *pceed)
403 {
404 DMCeed p = *pceed;
405
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 }
425
DMCeedComputeGeometry(DM dm,DMCeed sd)426 PetscErrorCode DMCeedComputeGeometry(DM dm, DMCeed sd)
427 {
428 #ifdef PETSC_HAVE_LIBCEED
429 Ceed ceed;
430 Vec coords;
431 CeedVector ccoords;
432 #endif
433
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 }
446
DMCeedComputeInfo(DM dm,DMCeed sd)447 PetscErrorCode DMCeedComputeInfo(DM dm, DMCeed sd)
448 {
449 #ifdef PETSC_HAVE_LIBCEED
450 CeedScalar *a;
451 #endif
452
453 PetscFunctionBegin;
454 #ifdef PETSC_HAVE_LIBCEED
455 PetscCallCEED(CeedVectorGetArrayWrite(sd->qi, CEED_MEM_HOST, &a));
456
457 IS iterIS;
458 DMLabel label = NULL;
459 const PetscInt *indices;
460 PetscInt value = 0, height = 1, NfInt = 0, Nf = 0;
461
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 }
473
474 PetscInt infoOffset = 0;
475
476 for (PetscInt p = 0; p < Nf; ++p) {
477 const PetscInt face = indices[p];
478 const PetscInt *supp;
479 PetscInt Ns;
480
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));
494
495 PetscCallCEED(CeedVectorRestoreArray(sd->qi, &a));
496 #endif
497 PetscFunctionReturn(PETSC_SUCCESS);
498 }
499