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(); 185 ierr = PetscOptionsGetEnum(NULL, NULL, "-dm_plex_shape", DMPlexShapes, (PetscEnum *) &shape, NULL);CHKERRQ(ierr); 186 switch (shape) { 187 case DM_SHAPE_BOX_SURFACE: 188 ctx->setupgeo = SetupMassGeoCube; 189 ctx->setupgeofname = SetupMassGeoCube_loc; 190 ctx->areaExact = 6.0; 191 break; 192 case DM_SHAPE_SPHERE: 193 ctx->setupgeo = SetupMassGeoSphere; 194 ctx->setupgeofname = SetupMassGeoSphere_loc; 195 ctx->areaExact = 4.0*M_PI; 196 break; 197 default: ctx->areaExact = 0.0; 198 } 199 200 ctx->apply = Mass; 201 ctx->applyfname = Mass_loc; 202 PetscFunctionReturn(0); 203 } 204 205 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) 206 { 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 211 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 212 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 213 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 214 PetscFunctionReturn(0); 215 } 216 217 static PetscErrorCode SetupDiscretization(DM dm) 218 { 219 DM cdm; 220 PetscFE fe; 221 PetscInt dim; 222 PetscBool simplex; 223 PetscErrorCode ierr; 224 225 PetscFunctionBeginUser; 226 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 227 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 228 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 229 ierr = PetscFESetName(fe, "indicator");CHKERRQ(ierr); 230 ierr = DMAddField(dm, NULL, (PetscObject) fe);CHKERRQ(ierr); 231 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 232 ierr = DMCreateDS(dm);CHKERRQ(ierr); 233 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL);CHKERRQ(ierr); 234 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 235 ierr = DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL);CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 239 static PetscErrorCode LibCeedSetupByDegree(DM dm, AppCtx *ctx, CeedData *data) 240 { 241 PetscDS ds; 242 PetscFE fe, cfe; 243 Ceed ceed; 244 CeedElemRestriction Erestrictx, Erestrictu, Erestrictq; 245 CeedQFunction qf_setupgeo; 246 CeedOperator op_setupgeo; 247 CeedVector xcoord; 248 CeedBasis basisu, basisx; 249 CeedInt Nqdata = 1; 250 CeedInt nqpts, nqptsx; 251 DM cdm; 252 Vec coords; 253 const PetscScalar *coordArray; 254 PetscInt dim, cdim, cStart, cEnd, Ncell; 255 PetscErrorCode ierr; 256 257 PetscFunctionBeginUser; 258 ierr = DMGetCeed(dm, &ceed);CHKERRQ(ierr); 259 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 260 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 261 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 262 Ncell = cEnd - cStart; 263 // CEED bases 264 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 265 ierr = PetscDSGetDiscretization(ds, 0, (PetscObject *) &fe);CHKERRQ(ierr); 266 ierr = PetscFEGetCeedBasis(fe, &basisu);CHKERRQ(ierr); 267 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 268 ierr = DMGetDS(cdm, &ds);CHKERRQ(ierr); 269 ierr = PetscDSGetDiscretization(ds, 0, (PetscObject *) &cfe);CHKERRQ(ierr); 270 ierr = PetscFEGetCeedBasis(cfe, &basisx);CHKERRQ(ierr); 271 272 ierr = DMPlexGetCeedRestriction(cdm, &Erestrictx);CHKERRQ(ierr); 273 ierr = DMPlexGetCeedRestriction(dm, &Erestrictu);CHKERRQ(ierr); 274 ierr = CeedBasisGetNumQuadraturePoints(basisu, &nqpts);CHKERRQ(ierr); 275 ierr = CeedBasisGetNumQuadraturePoints(basisx, &nqptsx);CHKERRQ(ierr); 276 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); 277 ierr = CeedElemRestrictionCreateStrided(ceed, Ncell, nqpts, Nqdata, Nqdata*Ncell*nqpts, CEED_STRIDES_BACKEND, &Erestrictq);CHKERRQ(ierr); 278 279 ierr = DMGetCoordinatesLocal(dm, &coords);CHKERRQ(ierr); 280 ierr = VecGetArrayRead(coords, &coordArray);CHKERRQ(ierr); 281 ierr = CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL);CHKERRQ(ierr); 282 ierr = CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *) coordArray);CHKERRQ(ierr); 283 ierr = VecRestoreArrayRead(coords, &coordArray);CHKERRQ(ierr); 284 285 // Create the vectors that will be needed in setup and apply 286 ierr = CeedElemRestrictionCreateVector(Erestrictu, &data->uceed, NULL);CHKERRQ(ierr); 287 ierr = CeedElemRestrictionCreateVector(Erestrictu, &data->vceed, NULL);CHKERRQ(ierr); 288 ierr = CeedElemRestrictionCreateVector(Erestrictq, &data->qdata, NULL);CHKERRQ(ierr); 289 290 // Create the Q-function that builds the operator (i.e. computes its quadrature data) and set its context data 291 ierr = CeedQFunctionCreateInterior(ceed, 1, ctx->setupgeo, ctx->setupgeofname, &qf_setupgeo);CHKERRQ(ierr); 292 ierr = CeedQFunctionAddInput(qf_setupgeo, "x", cdim, CEED_EVAL_INTERP);CHKERRQ(ierr); 293 ierr = CeedQFunctionAddInput(qf_setupgeo, "dx", cdim*dim, CEED_EVAL_GRAD);CHKERRQ(ierr); 294 ierr = CeedQFunctionAddInput(qf_setupgeo, "weight", 1, CEED_EVAL_WEIGHT);CHKERRQ(ierr); 295 ierr = CeedQFunctionAddOutput(qf_setupgeo, "qdata", Nqdata, CEED_EVAL_NONE);CHKERRQ(ierr); 296 297 // Set up the mass operator 298 ierr = CeedQFunctionCreateInterior(ceed, 1, ctx->apply, ctx->applyfname, &data->qf_apply);CHKERRQ(ierr); 299 ierr = CeedQFunctionAddInput(data->qf_apply, "u", 1, CEED_EVAL_INTERP);CHKERRQ(ierr); 300 ierr = CeedQFunctionAddInput(data->qf_apply, "qdata", Nqdata, CEED_EVAL_NONE);CHKERRQ(ierr); 301 ierr = CeedQFunctionAddOutput(data->qf_apply, "v", 1, CEED_EVAL_INTERP);CHKERRQ(ierr); 302 303 // Create the operator that builds the quadrature data for the operator 304 ierr = CeedOperatorCreate(ceed, qf_setupgeo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setupgeo);CHKERRQ(ierr); 305 ierr = CeedOperatorSetField(op_setupgeo, "x", Erestrictx, basisx, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 306 ierr = CeedOperatorSetField(op_setupgeo, "dx", Erestrictx, basisx, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 307 ierr = CeedOperatorSetField(op_setupgeo, "weight", CEED_ELEMRESTRICTION_NONE, basisx, CEED_VECTOR_NONE);CHKERRQ(ierr); 308 ierr = CeedOperatorSetField(op_setupgeo, "qdata", Erestrictq, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 309 310 // Create the mass operator 311 ierr = CeedOperatorCreate(ceed, data->qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &data->op_apply);CHKERRQ(ierr); 312 ierr = CeedOperatorSetField(data->op_apply, "u", Erestrictu, basisu, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 313 ierr = CeedOperatorSetField(data->op_apply, "qdata", Erestrictq, CEED_BASIS_COLLOCATED, data->qdata);CHKERRQ(ierr); 314 ierr = CeedOperatorSetField(data->op_apply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE);CHKERRQ(ierr); 315 316 // Setup qdata 317 ierr = CeedOperatorApply(op_setupgeo, xcoord, data->qdata, CEED_REQUEST_IMMEDIATE);CHKERRQ(ierr); 318 319 ierr = CeedElemRestrictionDestroy(&Erestrictq);CHKERRQ(ierr); 320 ierr = CeedQFunctionDestroy(&qf_setupgeo);CHKERRQ(ierr); 321 ierr = CeedOperatorDestroy(&op_setupgeo);CHKERRQ(ierr); 322 ierr = CeedVectorDestroy(&xcoord);CHKERRQ(ierr); 323 PetscFunctionReturn(0); 324 } 325 326 int main(int argc, char **argv) 327 { 328 MPI_Comm comm; 329 DM dm; 330 AppCtx ctx; 331 Vec U, Uloc, V, Vloc; 332 PetscScalar *v; 333 PetscScalar area; 334 CeedData ceeddata; 335 PetscErrorCode ierr; 336 337 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 338 comm = PETSC_COMM_WORLD; 339 ierr = ProcessOptions(comm, &ctx);CHKERRQ(ierr); 340 ierr = CreateMesh(comm, &ctx, &dm);CHKERRQ(ierr); 341 ierr = SetupDiscretization(dm);CHKERRQ(ierr); 342 343 ierr = LibCeedSetupByDegree(dm, &ctx, &ceeddata);CHKERRQ(ierr); 344 345 ierr = DMCreateGlobalVector(dm, &U);CHKERRQ(ierr); 346 ierr = DMCreateLocalVector(dm, &Uloc);CHKERRQ(ierr); 347 ierr = VecDuplicate(U, &V);CHKERRQ(ierr); 348 ierr = VecDuplicate(Uloc, &Vloc);CHKERRQ(ierr); 349 350 ierr = VecZeroEntries(V);CHKERRQ(ierr); 351 ierr = VecZeroEntries(Vloc);CHKERRQ(ierr); 352 ierr = VecGetArray(Vloc, &v);CHKERRQ(ierr); 353 ierr = CeedVectorSetArray(ceeddata.vceed, CEED_MEM_HOST, CEED_USE_POINTER, v);CHKERRQ(ierr); 354 ierr = CeedVectorSetValue(ceeddata.uceed, 1.0);CHKERRQ(ierr); 355 ierr = CeedOperatorApply(ceeddata.op_apply, ceeddata.uceed, ceeddata.vceed, CEED_REQUEST_IMMEDIATE);CHKERRQ(ierr); 356 ierr = CeedVectorTakeArray(ceeddata.vceed, CEED_MEM_HOST, NULL);CHKERRQ(ierr); 357 ierr = VecRestoreArray(Vloc, &v);CHKERRQ(ierr); 358 ierr = DMLocalToGlobalBegin(dm, Vloc, ADD_VALUES, V);CHKERRQ(ierr); 359 ierr = DMLocalToGlobalEnd(dm, Vloc, ADD_VALUES, V);CHKERRQ(ierr); 360 361 ierr = VecSum(V, &area);CHKERRQ(ierr); 362 if (ctx.areaExact > 0.) { 363 PetscReal error = PetscAbsReal(area - ctx.areaExact); 364 PetscReal tol = PETSC_SMALL; 365 366 ierr = PetscPrintf(comm, "Exact mesh surface area : % .14g\n", (double) ctx.areaExact);CHKERRQ(ierr); 367 ierr = PetscPrintf(comm, "Computed mesh surface area : % .14g\n", (double) area);CHKERRQ(ierr); 368 if (error > tol) { 369 ierr = PetscPrintf(comm, "Area error : % .14g\n", (double) error);CHKERRQ(ierr); 370 } else { 371 ierr = PetscPrintf(comm, "Area verifies!\n", (double) error);CHKERRQ(ierr); 372 } 373 } 374 375 ierr = CeedDataDestroy(&ceeddata);CHKERRQ(ierr); 376 ierr = VecDestroy(&U);CHKERRQ(ierr); 377 ierr = VecDestroy(&Uloc);CHKERRQ(ierr); 378 ierr = VecDestroy(&V);CHKERRQ(ierr); 379 ierr = VecDestroy(&Vloc);CHKERRQ(ierr); 380 ierr = DMDestroy(&dm);CHKERRQ(ierr); 381 return PetscFinalize(); 382 } 383 384 /*TEST 385 386 build: 387 requires: libceed 388 389 testset: 390 args: -dm_plex_simplex 0 -dm_distribute -petscspace_degree 3 -dm_view -dm_petscds_view \ 391 -petscfe_default_quadrature_order 4 -coord_dm_default_quadrature_order 4 392 393 test: 394 suffix: cube_3 395 args: -dm_plex_shape box_surface -dm_refine 2 396 397 test: 398 suffix: cube_3_p4 399 nsize: 4 400 args: -dm_refine_pre 1 -dm_plex_shape box_surface -dm_refine 1 401 402 test: 403 suffix: sphere_3 404 args: -dm_plex_shape sphere -dm_refine 3 405 406 test: 407 suffix: sphere_3_p4 408 nsize: 4 409 args: -dm_refine_pre 1 -dm_plex_shape sphere -dm_refine 2 410 411 TEST*/ 412