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*/