1 static char help[] = "Tests for geometry models\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petscds.h> 5 6 typedef struct { 7 PetscReal exactVol; // The exact volume of the shape 8 PetscReal volTol; // The relative tolerance for checking the volume 9 } AppCtx; 10 11 static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 12 { 13 f0[0] = 1.0; 14 } 15 16 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm) 17 { 18 PetscFunctionBegin; 19 PetscCall(DMCreate(comm, dm)); 20 PetscCall(DMSetType(*dm, DMPLEX)); 21 PetscCall(DMSetFromOptions(*dm)); 22 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 23 PetscFunctionReturn(PETSC_SUCCESS); 24 } 25 26 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 27 { 28 options->exactVol = 4. * PETSC_PI / 3.; 29 options->volTol = PETSC_SMALL; 30 31 PetscFunctionBeginUser; 32 PetscOptionsBegin(comm, "", "Geometry Model Test Options", "DMPLEX"); 33 PetscCall(PetscOptionsReal("-exact_vol", "Exact volume of the shape", __FILE__, options->exactVol, &options->exactVol, NULL)); 34 PetscCall(PetscOptionsReal("-vol_tol", "Relative tolerance for checking the volume", __FILE__, options->volTol, &options->volTol, NULL)); 35 PetscOptionsEnd(); 36 PetscFunctionReturn(PETSC_SUCCESS); 37 } 38 39 static PetscErrorCode CreateDiscretization(DM dm) 40 { 41 PetscFE fe; 42 PetscDS ds; 43 DMPolytopeType ct; 44 PetscInt dim, cStart; 45 46 PetscFunctionBeginUser; 47 PetscCall(DMGetDimension(dm, &dim)); 48 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 49 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 50 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, -1, &fe)); 51 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 52 PetscCall(DMCreateDS(dm)); 53 PetscCall(DMGetDS(dm, &ds)); 54 PetscCall(PetscDSSetObjective(ds, 0, identity)); 55 PetscCall(PetscFEDestroy(&fe)); 56 PetscFunctionReturn(PETSC_SUCCESS); 57 } 58 59 int main(int argc, char **argv) 60 { 61 DM dm; 62 Vec u; 63 PetscScalar volume; 64 AppCtx user; 65 66 PetscFunctionBeginUser; 67 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 68 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 69 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm)); 70 PetscCall(CreateDiscretization(dm)); 71 PetscCall(DMGetGlobalVector(dm, &u)); 72 PetscCall(VecSet(u, 0.)); 73 PetscCall(DMPlexComputeIntegralFEM(dm, u, &volume, NULL)); 74 PetscCall(DMRestoreGlobalVector(dm, &u)); 75 PetscCheck(PetscAbsScalar((volume - user.exactVol) / user.exactVol) < user.volTol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Invalid volume %g != %g", (double)PetscRealPart(volume), (double)user.exactVol); 76 PetscCall(DMDestroy(&dm)); 77 PetscCall(PetscFinalize()); 78 return 0; 79 } 80 81 /*TEST 82 83 # -dm_refine 6 -vol_tol 6e-4 works 84 test: 85 suffix: ball_0 86 requires: ctetgen 87 args: -dm_plex_dim 3 -dm_plex_shape ball -dm_refine 3 -dm_geom_model ball -vol_tol 5e-2 88 89 # -dm_refine 4 -vol_tol 2e-3 works 90 test: 91 suffix: cylinder_0 92 args: -dm_plex_dim 3 -dm_plex_shape cylinder -dm_plex_cylinder_num_refine 1 \ 93 -dm_refine 1 -dm_geom_model cylinder -exact_vol 3.141592653589 -vol_tol 1e-1 94 95 TEST*/ 96