1 static char help[] = "Tests 1D nested mesh refinement.\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petscds.h> 5 6 typedef struct { 7 PetscInt Nr; /* Number of refinements */ 8 PetscSimplePointFunc funcs[2]; /* Functions to test */ 9 } AppCtx; 10 11 static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 12 u[0] = 1.; 13 return 0; 14 } 15 16 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 17 u[0] = 2. * x[0] + 1.; 18 return 0; 19 } 20 21 static PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 22 u[0] = 3. * x[0] * x[0] + 2. * x[0] + 1.; 23 return 0; 24 } 25 26 static PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) { 27 u[0] = 4. * x[0] * x[0] * x[0] + 3. * x[0] * x[0] + 2. * x[0] + 1.; 28 return 0; 29 } 30 31 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 32 PetscFunctionBeginUser; 33 options->Nr = 1; 34 PetscOptionsBegin(comm, "", "1D Refinement Options", "DMPLEX"); 35 PetscCall(PetscOptionsInt("-num_refine", "Refine cycles", "ex46.c", options->Nr, &options->Nr, NULL)); 36 PetscOptionsEnd(); 37 PetscFunctionReturn(0); 38 } 39 40 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) { 41 PetscFunctionBeginUser; 42 PetscCall(DMCreate(comm, dm)); 43 PetscCall(DMSetType(*dm, DMPLEX)); 44 PetscCall(DMSetFromOptions(*dm)); 45 PetscCall(DMSetApplicationContext(*dm, user)); 46 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 47 PetscFunctionReturn(0); 48 } 49 50 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) { 51 DM cdm = dm; 52 PetscFE fe; 53 PetscSpace sp; 54 PetscInt dim, deg; 55 56 PetscFunctionBeginUser; 57 PetscCall(DMGetDimension(dm, &dim)); 58 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe)); 59 PetscCall(PetscObjectSetName((PetscObject)fe, "scalar")); 60 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 61 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe)); 62 PetscCall(DMCreateDS(dm)); 63 while (cdm) { 64 PetscCall(DMCopyDisc(dm, cdm)); 65 PetscCall(DMGetCoarseDM(cdm, &cdm)); 66 } 67 PetscCall(PetscFEGetBasisSpace(fe, &sp)); 68 PetscCall(PetscSpaceGetDegree(sp, °, NULL)); 69 switch (deg) { 70 case 0: user->funcs[0] = constant; break; 71 case 1: user->funcs[0] = linear; break; 72 case 2: user->funcs[0] = quadratic; break; 73 case 3: user->funcs[0] = cubic; break; 74 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine function to test for degree %" PetscInt_FMT, deg); 75 } 76 user->funcs[1] = user->funcs[0]; 77 PetscCall(PetscFEDestroy(&fe)); 78 PetscFunctionReturn(0); 79 } 80 81 static PetscErrorCode CheckError(DM dm, Vec u, PetscSimplePointFunc funcs[]) { 82 PetscReal error, tol = PETSC_SMALL; 83 MPI_Comm comm; 84 85 PetscFunctionBeginUser; 86 PetscCall(DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error)); 87 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 88 if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL at tolerance %g error %g\n", (double)tol, (double)error)); 89 else PetscCall(PetscPrintf(comm, "Function tests pass at tolerance %g\n", (double)tol)); 90 PetscFunctionReturn(0); 91 } 92 93 int main(int argc, char **argv) { 94 DM dm; 95 Vec u; 96 AppCtx user; 97 PetscInt cStart, cEnd, c, r; 98 99 PetscFunctionBeginUser; 100 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 101 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 102 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 103 PetscCall(SetupDiscretization(dm, &user)); 104 PetscCall(DMGetGlobalVector(dm, &u)); 105 PetscCall(DMProjectFunction(dm, 0.0, user.funcs, NULL, INSERT_ALL_VALUES, u)); 106 PetscCall(CheckError(dm, u, user.funcs)); 107 for (r = 0; r < user.Nr; ++r) { 108 DM adm; 109 DMLabel adapt; 110 Vec au; 111 Mat Interp; 112 113 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adapt)); 114 PetscCall(DMLabelSetDefaultValue(adapt, DM_ADAPT_COARSEN)); 115 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 116 for (c = cStart; c < cEnd; ++c) { 117 if (c % 2) PetscCall(DMLabelSetValue(adapt, c, DM_ADAPT_REFINE)); 118 } 119 PetscCall(DMAdaptLabel(dm, adapt, &adm)); 120 PetscCall(DMLabelDestroy(&adapt)); 121 PetscCall(PetscObjectSetName((PetscObject)adm, "Adapted Mesh")); 122 PetscCall(DMViewFromOptions(adm, NULL, "-dm_view")); 123 124 PetscCall(DMCreateInterpolation(dm, adm, &Interp, NULL)); 125 PetscCall(DMGetGlobalVector(adm, &au)); 126 PetscCall(MatInterpolate(Interp, u, au)); 127 PetscCall(CheckError(adm, au, user.funcs)); 128 PetscCall(MatDestroy(&Interp)); 129 PetscCall(DMRestoreGlobalVector(dm, &u)); 130 PetscCall(DMDestroy(&dm)); 131 dm = adm; 132 u = au; 133 } 134 PetscCall(DMRestoreGlobalVector(dm, &u)); 135 PetscCall(DMDestroy(&dm)); 136 PetscCall(PetscFinalize()); 137 return 0; 138 } 139 140 /*TEST 141 142 test: 143 suffix: 0 144 args: -num_refine 4 -petscspace_degree 3 \ 145 -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_plex_transform_type refine_1d -dm_plex_hash_location -dm_view 146 147 TEST*/ 148