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[1]; /* Functions to test */ 9 } AppCtx; 10 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 } 16 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 } 22 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 } 28 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 } 34 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(0); 43 } 44 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(0); 54 } 55 56 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 57 { 58 DM cdm = dm; 59 PetscFE fe; 60 PetscSpace sp; 61 PetscInt dim, deg; 62 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(DMCreateDS(dm)); 69 while (cdm) { 70 PetscCall(DMCopyDisc(dm,cdm)); 71 PetscCall(DMGetCoarseDM(cdm, &cdm)); 72 } 73 PetscCall(PetscFEGetBasisSpace(fe, &sp)); 74 PetscCall(PetscSpaceGetDegree(sp, °, NULL)); 75 switch (deg) { 76 case 0: user->funcs[0] = constant;break; 77 case 1: user->funcs[0] = linear;break; 78 case 2: user->funcs[0] = quadratic;break; 79 case 3: user->funcs[0] = cubic;break; 80 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Could not determine function to test for degree %" PetscInt_FMT, deg); 81 } 82 PetscCall(PetscFEDestroy(&fe)); 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode CheckError(DM dm, Vec u, PetscSimplePointFunc funcs[]) 87 { 88 PetscReal error, tol = PETSC_SMALL; 89 MPI_Comm comm; 90 91 PetscFunctionBeginUser; 92 PetscCall(DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error)); 93 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 94 if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL at tolerance %g error %g\n", (double)tol,(double) error)); 95 else PetscCall(PetscPrintf(comm, "Function tests pass at tolerance %g\n", (double)tol)); 96 PetscFunctionReturn(0); 97 } 98 99 int main(int argc, char **argv) 100 { 101 DM dm; 102 Vec u; 103 AppCtx user; 104 PetscInt cStart, cEnd, c, r; 105 106 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 107 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 108 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 109 PetscCall(SetupDiscretization(dm, &user)); 110 PetscCall(DMGetGlobalVector(dm, &u)); 111 PetscCall(DMProjectFunction(dm, 0.0, user.funcs, NULL, INSERT_ALL_VALUES, u)); 112 PetscCall(CheckError(dm, u, user.funcs)); 113 for (r = 0; r < user.Nr; ++r) { 114 DM adm; 115 DMLabel adapt; 116 Vec au; 117 Mat Interp; 118 119 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adapt)); 120 PetscCall(DMLabelSetDefaultValue(adapt, DM_ADAPT_COARSEN)); 121 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 122 for (c = cStart; c < cEnd; ++c) { 123 if (c % 2) PetscCall(DMLabelSetValue(adapt, c, DM_ADAPT_REFINE)); 124 } 125 PetscCall(DMAdaptLabel(dm, adapt, &adm)); 126 PetscCall(DMLabelDestroy(&adapt)); 127 PetscCall(PetscObjectSetName((PetscObject) adm, "Adapted Mesh")); 128 PetscCall(DMViewFromOptions(adm, NULL, "-dm_view")); 129 130 PetscCall(DMCreateInterpolation(dm, adm, &Interp, NULL)); 131 PetscCall(DMGetGlobalVector(adm, &au)); 132 PetscCall(MatInterpolate(Interp, u, au)); 133 PetscCall(CheckError(adm, au, user.funcs)); 134 PetscCall(MatDestroy(&Interp)); 135 PetscCall(DMRestoreGlobalVector(dm, &u)); 136 PetscCall(DMDestroy(&dm)); 137 dm = adm; 138 u = au; 139 } 140 PetscCall(DMRestoreGlobalVector(dm, &u)); 141 PetscCall(DMDestroy(&dm)); 142 PetscCall(PetscFinalize()); 143 return 0; 144 } 145 146 /*TEST 147 148 test: 149 suffix: 0 150 args: -num_refine 4 -petscspace_degree 3 \ 151 -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_plex_transform_type refine_1d -dm_plex_hash_location -dm_view 152 153 TEST*/ 154