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, cfe; 233 PetscInt dim, cnc; 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 = DMGetCoordinateDim(dm, &cnc);CHKERRQ(ierr); 247 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, cnc, simplex, NULL, PETSC_DETERMINE, &cfe);CHKERRQ(ierr); 248 ierr = DMProjectCoordinates(dm, cfe);CHKERRQ(ierr); 249 ierr = PetscFEDestroy(&cfe);CHKERRQ(ierr); 250 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 251 ierr = DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL);CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 static PetscErrorCode LibCeedSetupByDegree(DM dm, AppCtx *ctx, CeedData *data) 256 { 257 PetscDS ds; 258 PetscFE fe, cfe; 259 Ceed ceed; 260 CeedElemRestriction Erestrictx, Erestrictu, Erestrictq; 261 CeedQFunction qf_setupgeo; 262 CeedOperator op_setupgeo; 263 CeedVector xcoord; 264 CeedBasis basisu, basisx; 265 CeedInt Nqdata = 1; 266 CeedInt nqpts, nqptsx; 267 DM cdm; 268 Vec coords; 269 const PetscScalar *coordArray; 270 PetscInt dim, cdim, cStart, cEnd, Ncell; 271 PetscErrorCode ierr; 272 273 PetscFunctionBeginUser; 274 ierr = DMGetCeed(dm, &ceed);CHKERRQ(ierr); 275 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 276 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 277 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 278 Ncell = cEnd - cStart; 279 // CEED bases 280 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 281 ierr = PetscDSGetDiscretization(ds, 0, (PetscObject *) &fe);CHKERRQ(ierr); 282 ierr = PetscFEGetCeedBasis(fe, &basisu);CHKERRQ(ierr); 283 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 284 ierr = DMGetDS(cdm, &ds);CHKERRQ(ierr); 285 ierr = PetscDSGetDiscretization(ds, 0, (PetscObject *) &cfe);CHKERRQ(ierr); 286 ierr = PetscFEGetCeedBasis(cfe, &basisx);CHKERRQ(ierr); 287 288 ierr = DMPlexGetCeedRestriction(cdm, NULL, 0, 0, 0, &Erestrictx);CHKERRQ(ierr); 289 ierr = DMPlexGetCeedRestriction(dm, NULL, 0, 0, 0, &Erestrictu);CHKERRQ(ierr); 290 ierr = CeedBasisGetNumQuadraturePoints(basisu, &nqpts);CHKERRQ(ierr); 291 ierr = CeedBasisGetNumQuadraturePoints(basisx, &nqptsx);CHKERRQ(ierr); 292 PetscCheckFalse(nqptsx != nqpts,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of qpoints for u %D != %D Number of qpoints for x", nqpts, nqptsx); 293 ierr = CeedElemRestrictionCreateStrided(ceed, Ncell, nqpts, Nqdata, Nqdata*Ncell*nqpts, CEED_STRIDES_BACKEND, &Erestrictq);CHKERRQ(ierr); 294 295 ierr = DMGetCoordinatesLocal(dm, &coords);CHKERRQ(ierr); 296 ierr = VecGetArrayRead(coords, &coordArray);CHKERRQ(ierr); 297 ierr = CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL);CHKERRQ(ierr); 298 ierr = CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *) coordArray);CHKERRQ(ierr); 299 ierr = VecRestoreArrayRead(coords, &coordArray);CHKERRQ(ierr); 300 301 // Create the vectors that will be needed in setup and apply 302 ierr = CeedElemRestrictionCreateVector(Erestrictu, &data->uceed, NULL);CHKERRQ(ierr); 303 ierr = CeedElemRestrictionCreateVector(Erestrictu, &data->vceed, NULL);CHKERRQ(ierr); 304 ierr = CeedElemRestrictionCreateVector(Erestrictq, &data->qdata, NULL);CHKERRQ(ierr); 305 306 // Create the Q-function that builds the operator (i.e. computes its quadrature data) and set its context data 307 ierr = CeedQFunctionCreateInterior(ceed, 1, ctx->setupgeo, ctx->setupgeofname, &qf_setupgeo);CHKERRQ(ierr); 308 ierr = CeedQFunctionAddInput(qf_setupgeo, "x", cdim, CEED_EVAL_INTERP);CHKERRQ(ierr); 309 ierr = CeedQFunctionAddInput(qf_setupgeo, "dx", cdim*dim, CEED_EVAL_GRAD);CHKERRQ(ierr); 310 ierr = CeedQFunctionAddInput(qf_setupgeo, "weight", 1, CEED_EVAL_WEIGHT);CHKERRQ(ierr); 311 ierr = CeedQFunctionAddOutput(qf_setupgeo, "qdata", Nqdata, CEED_EVAL_NONE);CHKERRQ(ierr); 312 313 // Set up the mass operator 314 ierr = CeedQFunctionCreateInterior(ceed, 1, ctx->apply, ctx->applyfname, &data->qf_apply);CHKERRQ(ierr); 315 ierr = CeedQFunctionAddInput(data->qf_apply, "u", 1, CEED_EVAL_INTERP);CHKERRQ(ierr); 316 ierr = CeedQFunctionAddInput(data->qf_apply, "qdata", Nqdata, CEED_EVAL_NONE);CHKERRQ(ierr); 317 ierr = CeedQFunctionAddOutput(data->qf_apply, "v", 1, CEED_EVAL_INTERP);CHKERRQ(ierr); 318 319 // Create the operator that builds the quadrature data for the operator 320 ierr = CeedOperatorCreate(ceed, qf_setupgeo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setupgeo);CHKERRQ(ierr); 321 ierr = CeedOperatorSetField(op_setupgeo, "x", Erestrictx, basisx, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 322 ierr = CeedOperatorSetField(op_setupgeo, "dx", Erestrictx, basisx, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 323 ierr = CeedOperatorSetField(op_setupgeo, "weight", CEED_ELEMRESTRICTION_NONE, basisx, CEED_VECTOR_NONE);CHKERRQ(ierr); 324 ierr = CeedOperatorSetField(op_setupgeo, "qdata", Erestrictq, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 325 326 // Create the mass operator 327 ierr = CeedOperatorCreate(ceed, data->qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &data->op_apply);CHKERRQ(ierr); 328 ierr = CeedOperatorSetField(data->op_apply, "u", Erestrictu, basisu, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 329 ierr = CeedOperatorSetField(data->op_apply, "qdata", Erestrictq, CEED_BASIS_COLLOCATED, data->qdata);CHKERRQ(ierr); 330 ierr = CeedOperatorSetField(data->op_apply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 331 332 // Setup qdata 333 ierr = CeedOperatorApply(op_setupgeo, xcoord, data->qdata, CEED_REQUEST_IMMEDIATE);CHKERRQ(ierr); 334 335 ierr = CeedElemRestrictionDestroy(&Erestrictq);CHKERRQ(ierr); 336 ierr = CeedQFunctionDestroy(&qf_setupgeo);CHKERRQ(ierr); 337 ierr = CeedOperatorDestroy(&op_setupgeo);CHKERRQ(ierr); 338 ierr = CeedVectorDestroy(&xcoord);CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 int main(int argc, char **argv) 343 { 344 MPI_Comm comm; 345 DM dm; 346 AppCtx ctx; 347 Vec U, Uloc, V, Vloc; 348 PetscScalar *v; 349 PetscScalar area; 350 CeedData ceeddata; 351 PetscErrorCode ierr; 352 353 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 354 comm = PETSC_COMM_WORLD; 355 ierr = ProcessOptions(comm, &ctx);CHKERRQ(ierr); 356 ierr = CreateMesh(comm, &ctx, &dm);CHKERRQ(ierr); 357 ierr = SetupDiscretization(dm);CHKERRQ(ierr); 358 359 ierr = LibCeedSetupByDegree(dm, &ctx, &ceeddata);CHKERRQ(ierr); 360 361 ierr = DMCreateGlobalVector(dm, &U);CHKERRQ(ierr); 362 ierr = DMCreateLocalVector(dm, &Uloc);CHKERRQ(ierr); 363 ierr = VecDuplicate(U, &V);CHKERRQ(ierr); 364 ierr = VecDuplicate(Uloc, &Vloc);CHKERRQ(ierr); 365 366 /**/ 367 ierr = VecSet(Uloc, 1.);CHKERRQ(ierr); 368 ierr = VecZeroEntries(V);CHKERRQ(ierr); 369 ierr = VecZeroEntries(Vloc);CHKERRQ(ierr); 370 ierr = VecGetArray(Vloc, &v);CHKERRQ(ierr); 371 ierr = CeedVectorSetArray(ceeddata.vceed, CEED_MEM_HOST, CEED_USE_POINTER, v);CHKERRQ(ierr); 372 ierr = CeedVectorSetValue(ceeddata.uceed, 1.0);CHKERRQ(ierr); 373 ierr = CeedOperatorApply(ceeddata.op_apply, ceeddata.uceed, ceeddata.vceed, CEED_REQUEST_IMMEDIATE);CHKERRQ(ierr); 374 ierr = CeedVectorTakeArray(ceeddata.vceed, CEED_MEM_HOST, NULL);CHKERRQ(ierr); 375 ierr = VecRestoreArray(Vloc, &v);CHKERRQ(ierr); 376 ierr = DMLocalToGlobalBegin(dm, Vloc, ADD_VALUES, V);CHKERRQ(ierr); 377 ierr = DMLocalToGlobalEnd(dm, Vloc, ADD_VALUES, V);CHKERRQ(ierr); 378 379 ierr = VecSum(V, &area);CHKERRQ(ierr); 380 if (ctx.areaExact > 0.) { 381 PetscReal error = PetscAbsReal(area - ctx.areaExact); 382 PetscReal tol = PETSC_SMALL; 383 384 ierr = PetscPrintf(comm, "Exact mesh surface area : % .*f\n", fabs(ctx.areaExact - round(ctx.areaExact)) > 1E-15 ? 14 : 1, (double) ctx.areaExact);CHKERRQ(ierr); 385 ierr = PetscPrintf(comm, "Computed mesh surface area : % .*f\n", fabs(area - round(area)) > 1E-15 ? 14 : 1, (double) area);CHKERRQ(ierr); 386 if (error > tol) { 387 ierr = PetscPrintf(comm, "Area error : % .14f\n", (double) error);CHKERRQ(ierr); 388 } else { 389 ierr = PetscPrintf(comm, "Area verifies!\n", (double) error);CHKERRQ(ierr); 390 } 391 } 392 393 ierr = CeedDataDestroy(&ceeddata);CHKERRQ(ierr); 394 ierr = VecDestroy(&U);CHKERRQ(ierr); 395 ierr = VecDestroy(&Uloc);CHKERRQ(ierr); 396 ierr = VecDestroy(&V);CHKERRQ(ierr); 397 ierr = VecDestroy(&Vloc);CHKERRQ(ierr); 398 ierr = DMDestroy(&dm);CHKERRQ(ierr); 399 return PetscFinalize(); 400 } 401 402 /*TEST 403 404 build: 405 requires: libceed 406 407 testset: 408 args: -dm_plex_simplex 0 -petscspace_degree 3 -dm_view -dm_petscds_view \ 409 -petscfe_default_quadrature_order 4 -coord_dm_default_quadrature_order 4 410 411 test: 412 suffix: cube_3 413 args: -dm_plex_shape box_surface -dm_refine 2 414 415 test: 416 suffix: cube_3_p4 417 nsize: 4 418 args: -dm_refine_pre 1 -dm_plex_shape box_surface -dm_refine 1 419 420 test: 421 suffix: sphere_3 422 args: -dm_plex_shape sphere -dm_refine 3 423 424 test: 425 suffix: sphere_3_p4 426 nsize: 4 427 args: -dm_refine_pre 1 -dm_plex_shape sphere -dm_refine 2 428 429 TEST*/ 430