Actual source code: ex46.c

  1: static char help[] = "Tests 1D nested mesh refinement.\n\n";

  3: #include <petscdmplex.h>
  4: #include <petscds.h>

  6: typedef struct {
  7:   PetscInt             Nr;       /* Number of refinements */
  8:   PetscSimplePointFunc funcs[2]; /* Functions to test */
  9: } AppCtx;

 11: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 12: {
 13:   u[0] = 1.;
 14:   return PETSC_SUCCESS;
 15: }

 17: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 18: {
 19:   u[0] = 2. * x[0] + 1.;
 20:   return PETSC_SUCCESS;
 21: }

 23: static PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 24: {
 25:   u[0] = 3. * x[0] * x[0] + 2. * x[0] + 1.;
 26:   return PETSC_SUCCESS;
 27: }

 29: static PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 30: {
 31:   u[0] = 4. * x[0] * x[0] * x[0] + 3. * x[0] * x[0] + 2. * x[0] + 1.;
 32:   return PETSC_SUCCESS;
 33: }

 35: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 36: {
 37:   PetscFunctionBeginUser;
 38:   options->Nr = 1;
 39:   PetscOptionsBegin(comm, "", "1D Refinement Options", "DMPLEX");
 40:   PetscCall(PetscOptionsInt("-num_refine", "Refine cycles", "ex46.c", options->Nr, &options->Nr, NULL));
 41:   PetscOptionsEnd();
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

 45: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 46: {
 47:   PetscFunctionBeginUser;
 48:   PetscCall(DMCreate(comm, dm));
 49:   PetscCall(DMSetType(*dm, DMPLEX));
 50:   PetscCall(DMSetFromOptions(*dm));
 51:   PetscCall(DMSetApplicationContext(*dm, user));
 52:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
 57: {
 58:   DM         cdm = dm;
 59:   PetscFE    fe;
 60:   PetscSpace sp;
 61:   PetscInt   dim, deg;

 63:   PetscFunctionBeginUser;
 64:   PetscCall(DMGetDimension(dm, &dim));
 65:   PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe));
 66:   PetscCall(PetscObjectSetName((PetscObject)fe, "scalar"));
 67:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
 68:   PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe));
 69:   PetscCall(DMCreateDS(dm));
 70:   while (cdm) {
 71:     PetscCall(DMCopyDisc(dm, cdm));
 72:     PetscCall(DMGetCoarseDM(cdm, &cdm));
 73:   }
 74:   PetscCall(PetscFEGetBasisSpace(fe, &sp));
 75:   PetscCall(PetscSpaceGetDegree(sp, &deg, NULL));
 76:   switch (deg) {
 77:   case 0:
 78:     user->funcs[0] = constant;
 79:     break;
 80:   case 1:
 81:     user->funcs[0] = linear;
 82:     break;
 83:   case 2:
 84:     user->funcs[0] = quadratic;
 85:     break;
 86:   case 3:
 87:     user->funcs[0] = cubic;
 88:     break;
 89:   default:
 90:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine function to test for degree %" PetscInt_FMT, deg);
 91:   }
 92:   user->funcs[1] = user->funcs[0];
 93:   PetscCall(PetscFEDestroy(&fe));
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: static PetscErrorCode CheckError(DM dm, Vec u, PetscSimplePointFunc funcs[])
 98: {
 99:   PetscReal error, tol = PETSC_SMALL;
100:   MPI_Comm  comm;

102:   PetscFunctionBeginUser;
103:   PetscCall(DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error));
104:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
105:   if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL at tolerance %g error %g\n", (double)tol, (double)error));
106:   else PetscCall(PetscPrintf(comm, "Function tests pass at tolerance %g\n", (double)tol));
107:   PetscFunctionReturn(PETSC_SUCCESS);
108: }

110: int main(int argc, char **argv)
111: {
112:   DM       dm;
113:   Vec      u;
114:   AppCtx   user;
115:   PetscInt cStart, cEnd, c, r;

117:   PetscFunctionBeginUser;
118:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
119:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
120:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
121:   PetscCall(SetupDiscretization(dm, &user));
122:   PetscCall(DMGetGlobalVector(dm, &u));
123:   PetscCall(DMProjectFunction(dm, 0.0, user.funcs, NULL, INSERT_ALL_VALUES, u));
124:   PetscCall(CheckError(dm, u, user.funcs));
125:   for (r = 0; r < user.Nr; ++r) {
126:     DM      adm;
127:     DMLabel adapt;
128:     Vec     au;
129:     Mat     Interp;

131:     PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adapt));
132:     PetscCall(DMLabelSetDefaultValue(adapt, DM_ADAPT_COARSEN));
133:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
134:     for (c = cStart; c < cEnd; ++c) {
135:       if (c % 2) PetscCall(DMLabelSetValue(adapt, c, DM_ADAPT_REFINE));
136:     }
137:     PetscCall(DMAdaptLabel(dm, adapt, &adm));
138:     PetscCall(DMLabelDestroy(&adapt));
139:     PetscCall(PetscObjectSetName((PetscObject)adm, "Adapted Mesh"));
140:     PetscCall(DMViewFromOptions(adm, NULL, "-dm_view"));

142:     PetscCall(DMCreateInterpolation(dm, adm, &Interp, NULL));
143:     PetscCall(DMGetGlobalVector(adm, &au));
144:     PetscCall(MatInterpolate(Interp, u, au));
145:     PetscCall(CheckError(adm, au, user.funcs));
146:     PetscCall(MatDestroy(&Interp));
147:     PetscCall(DMRestoreGlobalVector(dm, &u));
148:     PetscCall(DMDestroy(&dm));
149:     dm = adm;
150:     u  = au;
151:   }
152:   PetscCall(DMRestoreGlobalVector(dm, &u));
153:   PetscCall(DMDestroy(&dm));
154:   PetscCall(PetscFinalize());
155:   return 0;
156: }

158: /*TEST

160:   test:
161:     suffix: 0
162:     args: -num_refine 4 -petscspace_degree 3 \
163:           -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_plex_transform_type refine_1d -dm_plex_hash_location -dm_view

165: TEST*/