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
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[])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
CreateMesh(MPI_Comm comm,DM * dm)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
ProcessOptions(MPI_Comm comm,AppCtx * options)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
CreateDiscretization(DM dm)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
main(int argc,char ** argv)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 output_file: output/empty.out
89
90 # -dm_refine 4 -vol_tol 2e-3 works
91 test:
92 suffix: cylinder_0
93 args: -dm_plex_dim 3 -dm_plex_shape cylinder -dm_plex_cylinder_num_refine 1 \
94 -dm_refine 1 -dm_geom_model cylinder -exact_vol 3.141592653589 -vol_tol 1e-1
95 output_file: output/empty.out
96
97 TEST*/
98