Actual source code: ex72.c
1: static char help[] = "Tests for geometry models\n\n";
3: #include <petscdmplex.h>
4: #include <petscds.h>
6: typedef struct {
7: PetscReal exactVol; // The exact volume of the shape
8: PetscReal volTol; // The relative tolerance for checking the volume
9: } AppCtx;
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: }
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: }
26: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
27: {
28: options->exactVol = 4. * PETSC_PI / 3.;
29: options->volTol = PETSC_SMALL;
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: }
39: static PetscErrorCode CreateDiscretization(DM dm)
40: {
41: PetscFE fe;
42: PetscDS ds;
43: DMPolytopeType ct;
44: PetscInt dim, cStart;
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: }
59: int main(int argc, char **argv)
60: {
61: DM dm;
62: Vec u;
63: PetscScalar volume;
64: AppCtx user;
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: }
81: /*TEST
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
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
95: TEST*/