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 0;
 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 0;
 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 0;
 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 0;
 33: }

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

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

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

 64:   DMGetDimension(dm, &dim);
 65:   PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe);
 66:   PetscObjectSetName((PetscObject)fe, "scalar");
 67:   DMSetField(dm, 0, NULL, (PetscObject)fe);
 68:   DMSetField(dm, 1, NULL, (PetscObject)fe);
 69:   DMCreateDS(dm);
 70:   while (cdm) {
 71:     DMCopyDisc(dm, cdm);
 72:     DMGetCoarseDM(cdm, &cdm);
 73:   }
 74:   PetscFEGetBasisSpace(fe, &sp);
 75:   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:   PetscFEDestroy(&fe);
 94:   return 0;
 95: }

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

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

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

118:   PetscInitialize(&argc, &argv, NULL, help);
119:   ProcessOptions(PETSC_COMM_WORLD, &user);
120:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
121:   SetupDiscretization(dm, &user);
122:   DMGetGlobalVector(dm, &u);
123:   DMProjectFunction(dm, 0.0, user.funcs, NULL, INSERT_ALL_VALUES, u);
124:   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:     DMLabelCreate(PETSC_COMM_SELF, "adapt", &adapt);
132:     DMLabelSetDefaultValue(adapt, DM_ADAPT_COARSEN);
133:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
134:     for (c = cStart; c < cEnd; ++c) {
135:       if (c % 2) DMLabelSetValue(adapt, c, DM_ADAPT_REFINE);
136:     }
137:     DMAdaptLabel(dm, adapt, &adm);
138:     DMLabelDestroy(&adapt);
139:     PetscObjectSetName((PetscObject)adm, "Adapted Mesh");
140:     DMViewFromOptions(adm, NULL, "-dm_view");

142:     DMCreateInterpolation(dm, adm, &Interp, NULL);
143:     DMGetGlobalVector(adm, &au);
144:     MatInterpolate(Interp, u, au);
145:     CheckError(adm, au, user.funcs);
146:     MatDestroy(&Interp);
147:     DMRestoreGlobalVector(dm, &u);
148:     DMDestroy(&dm);
149:     dm = adm;
150:     u  = au;
151:   }
152:   DMRestoreGlobalVector(dm, &u);
153:   DMDestroy(&dm);
154:   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*/