xref: /petsc/src/dm/impls/plex/tests/ex42.c (revision dfd57a172ac9fa6c7b5fe6de6ab5df85cefc2996)
1 static const char help[] = "Simple libCEED test to calculate surface area using 1^T M 1";
2 
3 /*
4   This is a recreation of libCeed Example 2: https://libceed.readthedocs.io/en/latest/examples/ceed/
5 */
6 
7 #include <petscdmceed.h>
8 #include <petscdmplexceed.h>
9 #include <petscfeceed.h>
10 #include <petscdmplex.h>
11 #include <petscds.h>
12 
13 typedef struct {
14   PetscReal         areaExact;
15   CeedQFunctionUser setupgeo,       apply;
16   const char       *setupgeofname, *applyfname;
17 } AppCtx;
18 
19 typedef struct {
20   CeedQFunction qf_apply;
21   CeedOperator  op_apply;
22   CeedVector    qdata, uceed, vceed;
23 } CeedData;
24 
25 static PetscErrorCode CeedDataDestroy(CeedData *data)
26 {
27   PetscErrorCode ierr;
28 
29   PetscFunctionBeginUser;
30   ierr = CeedVectorDestroy(&data->qdata);CHKERRQ(ierr);
31   ierr = CeedVectorDestroy(&data->uceed);CHKERRQ(ierr);
32   ierr = CeedVectorDestroy(&data->vceed);CHKERRQ(ierr);
33   ierr = CeedQFunctionDestroy(&data->qf_apply);CHKERRQ(ierr);
34   ierr = CeedOperatorDestroy(&data->op_apply);CHKERRQ(ierr);
35   PetscFunctionReturn(0);
36 }
37 
38 CEED_QFUNCTION(Mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
39 {
40   const CeedScalar *u = in[0], *qdata = in[1];
41   CeedScalar       *v = out[0];
42 
43   CeedPragmaSIMD
44   for (CeedInt i = 0; i < Q; ++i)
45     v[i] = qdata[i] * u[i];
46 
47   return 0;
48 }
49 
50 /*
51 // Reference (parent) 2D coordinates: X \in [-1, 1]^2
52 //
53 // Global physical coordinates given by the mesh (3D): xx \in [-l, l]^3
54 //
55 // Local physical coordinates on the manifold (2D): x \in [-l, l]^2
56 //
57 // Change of coordinates matrix computed by the library:
58 //   (physical 3D coords relative to reference 2D coords)
59 //   dxx_j/dX_i (indicial notation) [3 * 2]
60 //
61 // Change of coordinates x (physical 2D) relative to xx (phyisical 3D):
62 //   dx_i/dxx_j (indicial notation) [2 * 3]
63 //
64 // Change of coordinates x (physical 2D) relative to X (reference 2D):
65 //   (by chain rule)
66 //   dx_i/dX_j = dx_i/dxx_k * dxx_k/dX_j
67 //
68 // The quadrature data is stored in the array qdata.
69 //
70 // We require the determinant of the Jacobian to properly compute integrals of the form: int(u v)
71 //
72 // Qdata: w * det(dx_i/dX_j)
73 */
74 CEED_QFUNCTION(SetupMassGeoCube)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
75 {
76   const CeedScalar *J = in[1], *w = in[2];
77   CeedScalar       *qdata = out[0];
78 
79   CeedPragmaSIMD
80   for (CeedInt i = 0; i < Q; ++i) {
81     // Read dxxdX Jacobian entries, stored as [[0 3], [1 4], [2 5]]
82     const CeedScalar dxxdX[3][2] = {{J[i+Q*0], J[i+Q*3]},
83                                     {J[i+Q*1], J[i+Q*4]},
84                                     {J[i+Q*2], J[i+Q*5]}};
85     // Modulus of dxxdX column vectors
86     const CeedScalar modg1 = PetscSqrtReal(dxxdX[0][0]*dxxdX[0][0] + dxxdX[1][0]*dxxdX[1][0] + dxxdX[2][0]*dxxdX[2][0]);
87     const CeedScalar modg2 = PetscSqrtReal(dxxdX[0][1]*dxxdX[0][1] + dxxdX[1][1]*dxxdX[1][1] + dxxdX[2][1]*dxxdX[2][1]);
88     // Use normalized column vectors of dxxdX as rows of dxdxx
89     const CeedScalar dxdxx[2][3] = {{dxxdX[0][0] / modg1, dxxdX[1][0] / modg1, dxxdX[2][0] / modg1},
90                                     {dxxdX[0][1] / modg2, dxxdX[1][1] / modg2, dxxdX[2][1] / modg2}};
91 
92     CeedScalar dxdX[2][2];
93     for (int j = 0; j < 2; ++j)
94       for (int k = 0; k < 2; ++k) {
95         dxdX[j][k] = 0;
96         for (int l = 0; l < 3; ++l)
97           dxdX[j][k] += dxdxx[j][l]*dxxdX[l][k];
98       }
99     qdata[i+Q*0] = (dxdX[0][0]*dxdX[1][1] - dxdX[1][0]*dxdX[0][1]) * w[i]; /* det J * weight */
100   }
101   return 0;
102 }
103 
104 /*
105 // Reference (parent) 2D coordinates: X \in [-1, 1]^2
106 //
107 // Global 3D physical coordinates given by the mesh: xx \in [-R, R]^3
108 //   with R radius of the sphere
109 //
110 // Local 3D physical coordinates on the 2D manifold: x \in [-l, l]^3
111 //   with l half edge of the cube inscribed in the sphere
112 //
113 // Change of coordinates matrix computed by the library:
114 //   (physical 3D coords relative to reference 2D coords)
115 //   dxx_j/dX_i (indicial notation) [3 * 2]
116 //
117 // Change of coordinates x (on the 2D manifold) relative to xx (phyisical 3D):
118 //   dx_i/dxx_j (indicial notation) [3 * 3]
119 //
120 // Change of coordinates x (on the 2D manifold) relative to X (reference 2D):
121 //   (by chain rule)
122 //   dx_i/dX_j = dx_i/dxx_k * dxx_k/dX_j [3 * 2]
123 //
124 // modJ is given by the magnitude of the cross product of the columns of dx_i/dX_j
125 //
126 // The quadrature data is stored in the array qdata.
127 //
128 // We require the determinant of the Jacobian to properly compute integrals of
129 //   the form: int(u v)
130 //
131 // Qdata: modJ * w
132 */
133 CEED_QFUNCTION(SetupMassGeoSphere)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
134   const CeedScalar *X = in[0], *J = in[1], *w = in[2];
135   CeedScalar       *qdata = out[0];
136 
137   CeedPragmaSIMD
138   for (CeedInt i = 0; i < Q; ++i) {
139     const CeedScalar xx[3][1] = {{X[i+0*Q]}, {X[i+1*Q]}, {X[i+2*Q]}};
140     // Read dxxdX Jacobian entries, stored as [[0 3], [1 4], [2 5]]
141     const CeedScalar dxxdX[3][2] = {{J[i+Q*0], J[i+Q*3]},
142                                     {J[i+Q*1], J[i+Q*4]},
143                                     {J[i+Q*2], J[i+Q*5]}};
144     // Setup
145     const CeedScalar modxxsq = xx[0][0]*xx[0][0]+xx[1][0]*xx[1][0]+xx[2][0]*xx[2][0];
146     CeedScalar xxsq[3][3];
147     for (int j = 0; j < 3; ++j)
148       for (int k = 0; k < 3; ++k) {
149         xxsq[j][k] = 0.;
150         for (int l = 0; l < 1; ++l)
151           xxsq[j][k] += xx[j][l]*xx[k][l] / (sqrt(modxxsq) * modxxsq);
152       }
153 
154     const CeedScalar dxdxx[3][3] = {{1./sqrt(modxxsq) - xxsq[0][0], -xxsq[0][1], -xxsq[0][2]},
155                                     {-xxsq[1][0], 1./sqrt(modxxsq) - xxsq[1][1], -xxsq[1][2]},
156                                     {-xxsq[2][0], -xxsq[2][1], 1./sqrt(modxxsq) - xxsq[2][2]}};
157 
158     CeedScalar dxdX[3][2];
159     for (int j = 0; j < 3; ++j)
160       for (int k = 0; k < 2; ++k) {
161         dxdX[j][k] = 0.;
162         for (int l = 0; l < 3; ++l)
163           dxdX[j][k] += dxdxx[j][l]*dxxdX[l][k];
164       }
165     // J is given by the cross product of the columns of dxdX
166     const CeedScalar J[3][1] = {{dxdX[1][0]*dxdX[2][1] - dxdX[2][0]*dxdX[1][1]},
167                                 {dxdX[2][0]*dxdX[0][1] - dxdX[0][0]*dxdX[2][1]},
168                                 {dxdX[0][0]*dxdX[1][1] - dxdX[1][0]*dxdX[0][1]}};
169     // Use the magnitude of J as our detJ (volume scaling factor)
170     const CeedScalar modJ = sqrt(J[0][0]*J[0][0]+J[1][0]*J[1][0]+J[2][0]*J[2][0]);
171     qdata[i+Q*0] = modJ * w[i];
172   }
173   return 0;
174 }
175 
176 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *ctx)
177 {
178   DMPlexShape    shape = DM_SHAPE_UNKNOWN;
179   PetscErrorCode ierr;
180 
181   PetscFunctionBeginUser;
182 
183   ierr = PetscOptionsBegin(comm, "", "libCEED Test Options", "DMPLEX");CHKERRQ(ierr);
184   ierr = PetscOptionsEnd();CHKERRQ(ierr);
185   ierr = PetscOptionsGetEnum(NULL, NULL, "-dm_plex_shape", DMPlexShapes, (PetscEnum *) &shape, NULL);CHKERRQ(ierr);
186   ctx->setupgeo      = NULL;
187   ctx->setupgeofname = NULL;
188   ctx->apply         = Mass;
189   ctx->applyfname    = Mass_loc;
190   ctx->areaExact     = 0.0;
191   switch (shape) {
192     case DM_SHAPE_BOX_SURFACE:
193       ctx->setupgeo      = SetupMassGeoCube;
194       ctx->setupgeofname = SetupMassGeoCube_loc;
195       ctx->areaExact     = 6.0;
196       break;
197     case DM_SHAPE_SPHERE:
198       ctx->setupgeo      = SetupMassGeoSphere;
199       ctx->setupgeofname = SetupMassGeoSphere_loc;
200       ctx->areaExact     = 4.0*M_PI;
201       break;
202     default: break;
203   }
204   PetscFunctionReturn(0);
205 }
206 
207 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
208 {
209   PetscErrorCode ierr;
210 
211   PetscFunctionBegin;
212   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
213   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
214   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
215   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
216 #ifdef PETSC_HAVE_LIBCEED
217   {
218     Ceed        ceed;
219     const char *usedresource;
220 
221     ierr = DMGetCeed(*dm, &ceed);CHKERRQ(ierr);
222     ierr = CeedGetResource(ceed, &usedresource);CHKERRQ(ierr);
223     ierr = PetscPrintf(PetscObjectComm((PetscObject) *dm), "libCEED Backend: %s\n", usedresource);CHKERRQ(ierr);
224   }
225 #endif
226   PetscFunctionReturn(0);
227 }
228 
229 static PetscErrorCode SetupDiscretization(DM dm)
230 {
231   DM             cdm;
232   PetscFE        fe;
233   PetscInt       dim;
234   PetscBool      simplex;
235   PetscErrorCode ierr;
236 
237   PetscFunctionBeginUser;
238   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
239   ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr);
240   ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fe);CHKERRQ(ierr);
241   ierr = PetscFESetName(fe, "indicator");CHKERRQ(ierr);
242   ierr = DMAddField(dm, NULL, (PetscObject) fe);CHKERRQ(ierr);
243   ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
244   ierr = DMCreateDS(dm);CHKERRQ(ierr);
245   ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL);CHKERRQ(ierr);
246   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
247   ierr = DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL);CHKERRQ(ierr);
248   PetscFunctionReturn(0);
249 }
250 
251 static PetscErrorCode LibCeedSetupByDegree(DM dm, AppCtx *ctx, CeedData *data)
252 {
253   PetscDS             ds;
254   PetscFE             fe, cfe;
255   Ceed                ceed;
256   CeedElemRestriction Erestrictx, Erestrictu, Erestrictq;
257   CeedQFunction       qf_setupgeo;
258   CeedOperator        op_setupgeo;
259   CeedVector          xcoord;
260   CeedBasis           basisu, basisx;
261   CeedInt             Nqdata = 1;
262   CeedInt             nqpts, nqptsx;
263   DM                  cdm;
264   Vec                 coords;
265   const PetscScalar  *coordArray;
266   PetscInt            dim, cdim, cStart, cEnd, Ncell;
267   PetscErrorCode      ierr;
268 
269   PetscFunctionBeginUser;
270   ierr = DMGetCeed(dm, &ceed);CHKERRQ(ierr);
271   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
272   ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr);
273   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
274   Ncell = cEnd - cStart;
275   // CEED bases
276   ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
277   ierr = PetscDSGetDiscretization(ds, 0, (PetscObject *) &fe);CHKERRQ(ierr);
278   ierr = PetscFEGetCeedBasis(fe, &basisu);CHKERRQ(ierr);
279   ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
280   ierr = DMGetDS(cdm, &ds);CHKERRQ(ierr);
281   ierr = PetscDSGetDiscretization(ds, 0, (PetscObject *) &cfe);CHKERRQ(ierr);
282   ierr = PetscFEGetCeedBasis(cfe, &basisx);CHKERRQ(ierr);
283 
284   ierr = DMPlexGetCeedRestriction(cdm, &Erestrictx);CHKERRQ(ierr);
285   ierr = DMPlexGetCeedRestriction(dm,  &Erestrictu);CHKERRQ(ierr);
286   ierr = CeedBasisGetNumQuadraturePoints(basisu, &nqpts);CHKERRQ(ierr);
287   ierr = CeedBasisGetNumQuadraturePoints(basisx, &nqptsx);CHKERRQ(ierr);
288   if (nqptsx != nqpts) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of qpoints for u %D != %D Number of qpoints for x", nqpts, nqptsx);
289   ierr = CeedElemRestrictionCreateStrided(ceed, Ncell, nqpts, Nqdata, Nqdata*Ncell*nqpts, CEED_STRIDES_BACKEND, &Erestrictq);CHKERRQ(ierr);
290 
291   ierr = DMGetCoordinatesLocal(dm, &coords);CHKERRQ(ierr);
292   ierr = VecGetArrayRead(coords, &coordArray);CHKERRQ(ierr);
293   ierr = CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL);CHKERRQ(ierr);
294   ierr = CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *) coordArray);CHKERRQ(ierr);
295   ierr = VecRestoreArrayRead(coords, &coordArray);CHKERRQ(ierr);
296 
297   // Create the vectors that will be needed in setup and apply
298   ierr = CeedElemRestrictionCreateVector(Erestrictu, &data->uceed, NULL);CHKERRQ(ierr);
299   ierr = CeedElemRestrictionCreateVector(Erestrictu, &data->vceed, NULL);CHKERRQ(ierr);
300   ierr = CeedElemRestrictionCreateVector(Erestrictq, &data->qdata, NULL);CHKERRQ(ierr);
301 
302   // Create the Q-function that builds the operator (i.e. computes its quadrature data) and set its context data
303   ierr = CeedQFunctionCreateInterior(ceed, 1, ctx->setupgeo, ctx->setupgeofname, &qf_setupgeo);CHKERRQ(ierr);
304   ierr = CeedQFunctionAddInput(qf_setupgeo,  "x",      cdim,     CEED_EVAL_INTERP);CHKERRQ(ierr);
305   ierr = CeedQFunctionAddInput(qf_setupgeo,  "dx",     cdim*dim, CEED_EVAL_GRAD);CHKERRQ(ierr);
306   ierr = CeedQFunctionAddInput(qf_setupgeo,  "weight", 1,        CEED_EVAL_WEIGHT);CHKERRQ(ierr);
307   ierr = CeedQFunctionAddOutput(qf_setupgeo, "qdata",  Nqdata,   CEED_EVAL_NONE);CHKERRQ(ierr);
308 
309   // Set up the mass operator
310   ierr = CeedQFunctionCreateInterior(ceed, 1, ctx->apply, ctx->applyfname, &data->qf_apply);CHKERRQ(ierr);
311   ierr = CeedQFunctionAddInput(data->qf_apply,  "u",     1,      CEED_EVAL_INTERP);CHKERRQ(ierr);
312   ierr = CeedQFunctionAddInput(data->qf_apply,  "qdata", Nqdata, CEED_EVAL_NONE);CHKERRQ(ierr);
313   ierr = CeedQFunctionAddOutput(data->qf_apply, "v",     1,      CEED_EVAL_INTERP);CHKERRQ(ierr);
314 
315   // Create the operator that builds the quadrature data for the operator
316   ierr = CeedOperatorCreate(ceed, qf_setupgeo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setupgeo);CHKERRQ(ierr);
317   ierr = CeedOperatorSetField(op_setupgeo, "x",      Erestrictx, basisx, CEED_VECTOR_ACTIVE);CHKERRQ(ierr);
318   ierr = CeedOperatorSetField(op_setupgeo, "dx",     Erestrictx, basisx, CEED_VECTOR_ACTIVE);CHKERRQ(ierr);
319   ierr = CeedOperatorSetField(op_setupgeo, "weight", CEED_ELEMRESTRICTION_NONE, basisx, CEED_VECTOR_NONE);CHKERRQ(ierr);
320   ierr = CeedOperatorSetField(op_setupgeo, "qdata",  Erestrictq, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);CHKERRQ(ierr);
321 
322   // Create the mass operator
323   ierr = CeedOperatorCreate(ceed, data->qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &data->op_apply);CHKERRQ(ierr);
324   ierr = CeedOperatorSetField(data->op_apply, "u",     Erestrictu, basisu, CEED_VECTOR_ACTIVE);CHKERRQ(ierr);
325   ierr = CeedOperatorSetField(data->op_apply, "qdata", Erestrictq, CEED_BASIS_COLLOCATED, data->qdata);CHKERRQ(ierr);
326   ierr = CeedOperatorSetField(data->op_apply, "v",     Erestrictu, basisu, CEED_VECTOR_ACTIVE);CHKERRQ(ierr);
327 
328   // Setup qdata
329   ierr = CeedOperatorApply(op_setupgeo, xcoord, data->qdata, CEED_REQUEST_IMMEDIATE);CHKERRQ(ierr);
330 
331   ierr = CeedElemRestrictionDestroy(&Erestrictq);CHKERRQ(ierr);
332   ierr = CeedQFunctionDestroy(&qf_setupgeo);CHKERRQ(ierr);
333   ierr = CeedOperatorDestroy(&op_setupgeo);CHKERRQ(ierr);
334   ierr = CeedVectorDestroy(&xcoord);CHKERRQ(ierr);
335   PetscFunctionReturn(0);
336 }
337 
338 int main(int argc, char **argv)
339 {
340   MPI_Comm       comm;
341   DM             dm;
342   AppCtx         ctx;
343   Vec            U, Uloc, V, Vloc;
344   PetscScalar   *v;
345   PetscScalar    area;
346   CeedData       ceeddata;
347   PetscErrorCode ierr;
348 
349   ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
350   comm = PETSC_COMM_WORLD;
351   ierr = ProcessOptions(comm, &ctx);CHKERRQ(ierr);
352   ierr = CreateMesh(comm, &ctx, &dm);CHKERRQ(ierr);
353   ierr = SetupDiscretization(dm);CHKERRQ(ierr);
354 
355   ierr = LibCeedSetupByDegree(dm, &ctx, &ceeddata);CHKERRQ(ierr);
356 
357   ierr = DMCreateGlobalVector(dm, &U);CHKERRQ(ierr);
358   ierr = DMCreateLocalVector(dm, &Uloc);CHKERRQ(ierr);
359   ierr = VecDuplicate(U, &V);CHKERRQ(ierr);
360   ierr = VecDuplicate(Uloc, &Vloc);CHKERRQ(ierr);
361 
362   /**/
363   ierr = VecZeroEntries(V);CHKERRQ(ierr);
364   ierr = VecZeroEntries(Vloc);CHKERRQ(ierr);
365   ierr = VecGetArray(Vloc, &v);CHKERRQ(ierr);
366   ierr = CeedVectorSetArray(ceeddata.vceed, CEED_MEM_HOST, CEED_USE_POINTER, v);CHKERRQ(ierr);
367   ierr = CeedVectorSetValue(ceeddata.uceed, 1.0);CHKERRQ(ierr);
368   ierr = CeedOperatorApply(ceeddata.op_apply, ceeddata.uceed, ceeddata.vceed, CEED_REQUEST_IMMEDIATE);CHKERRQ(ierr);
369   ierr = CeedVectorTakeArray(ceeddata.vceed, CEED_MEM_HOST, NULL);CHKERRQ(ierr);
370   ierr = VecRestoreArray(Vloc, &v);CHKERRQ(ierr);
371   ierr = DMLocalToGlobalBegin(dm, Vloc, ADD_VALUES, V);CHKERRQ(ierr);
372   ierr = DMLocalToGlobalEnd(dm, Vloc, ADD_VALUES, V);CHKERRQ(ierr);
373 
374   ierr = VecSum(V, &area);CHKERRQ(ierr);
375   if (ctx.areaExact > 0.) {
376     PetscReal error = PetscAbsReal(area - ctx.areaExact);
377     PetscReal tol   = PETSC_SMALL;
378 
379     ierr = PetscPrintf(comm,   "Exact mesh surface area    : % .14g\n", (double) ctx.areaExact);CHKERRQ(ierr);
380     ierr = PetscPrintf(comm,   "Computed mesh surface area : % .14g\n", (double) area);CHKERRQ(ierr);
381     if (error > tol) {
382       ierr = PetscPrintf(comm, "Area error                 : % .14g\n", (double) error);CHKERRQ(ierr);
383     } else {
384       ierr = PetscPrintf(comm, "Area verifies!\n", (double) error);CHKERRQ(ierr);
385     }
386   }
387 
388   ierr = CeedDataDestroy(&ceeddata);CHKERRQ(ierr);
389   ierr = VecDestroy(&U);CHKERRQ(ierr);
390   ierr = VecDestroy(&Uloc);CHKERRQ(ierr);
391   ierr = VecDestroy(&V);CHKERRQ(ierr);
392   ierr = VecDestroy(&Vloc);CHKERRQ(ierr);
393   ierr = DMDestroy(&dm);CHKERRQ(ierr);
394   return PetscFinalize();
395 }
396 
397 /*TEST
398 
399   build:
400     requires: libceed
401 
402   testset:
403     args: -dm_plex_simplex 0 -dm_distribute -petscspace_degree 3 -dm_view -dm_petscds_view \
404           -petscfe_default_quadrature_order 4 -coord_dm_default_quadrature_order 4
405 
406     test:
407       suffix: cube_3
408       args: -dm_plex_shape box_surface -dm_refine 2
409 
410     test:
411       suffix: cube_3_p4
412       nsize: 4
413       args: -dm_refine_pre 1 -dm_plex_shape box_surface -dm_refine 1
414 
415     test:
416       suffix: sphere_3
417       args: -dm_plex_shape sphere -dm_refine 3
418 
419     test:
420       suffix: sphere_3_p4
421       nsize: 4
422       args: -dm_refine_pre 1 -dm_plex_shape sphere -dm_refine 2
423 
424 TEST*/
425