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