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");CHKERRQ(ierr); 42 CHKERRQ(PetscOptionsInt("-num_refine", "Refine cycles", "ex46.c", options->Nr, &options->Nr, NULL)); 43 ierr = PetscOptionsEnd();CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 48 { 49 PetscFunctionBeginUser; 50 CHKERRQ(DMCreate(comm, dm)); 51 CHKERRQ(DMSetType(*dm, DMPLEX)); 52 CHKERRQ(DMSetFromOptions(*dm)); 53 CHKERRQ(DMSetApplicationContext(*dm, user)); 54 CHKERRQ(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 CHKERRQ(DMGetDimension(dm, &dim)); 67 CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe)); 68 CHKERRQ(PetscObjectSetName((PetscObject) fe, "scalar")); 69 CHKERRQ(DMSetField(dm, 0, NULL, (PetscObject) fe)); 70 CHKERRQ(DMCreateDS(dm)); 71 while (cdm) { 72 CHKERRQ(DMCopyDisc(dm,cdm)); 73 CHKERRQ(DMGetCoarseDM(cdm, &cdm)); 74 } 75 CHKERRQ(PetscFEGetBasisSpace(fe, &sp)); 76 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error)); 95 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 96 if (error > tol) CHKERRQ(PetscPrintf(comm, "Function tests FAIL at tolerance %g error %g\n", (double)tol,(double) error)); 97 else CHKERRQ(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 PetscErrorCode ierr; 108 109 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 110 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 111 CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 112 CHKERRQ(SetupDiscretization(dm, &user)); 113 CHKERRQ(DMGetGlobalVector(dm, &u)); 114 CHKERRQ(DMProjectFunction(dm, 0.0, user.funcs, NULL, INSERT_ALL_VALUES, u)); 115 CHKERRQ(CheckError(dm, u, user.funcs)); 116 for (r = 0; r < user.Nr; ++r) { 117 DM adm; 118 DMLabel adapt; 119 Vec au; 120 Mat Interp; 121 122 CHKERRQ(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adapt)); 123 CHKERRQ(DMLabelSetDefaultValue(adapt, DM_ADAPT_COARSEN)); 124 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 125 for (c = cStart; c < cEnd; ++c) { 126 if (c % 2) CHKERRQ(DMLabelSetValue(adapt, c, DM_ADAPT_REFINE)); 127 } 128 CHKERRQ(DMAdaptLabel(dm, adapt, &adm)); 129 CHKERRQ(DMLabelDestroy(&adapt)); 130 CHKERRQ(PetscObjectSetName((PetscObject) adm, "Adapted Mesh")); 131 CHKERRQ(DMViewFromOptions(adm, NULL, "-dm_view")); 132 133 CHKERRQ(DMCreateInterpolation(dm, adm, &Interp, NULL)); 134 CHKERRQ(DMGetGlobalVector(adm, &au)); 135 CHKERRQ(MatInterpolate(Interp, u, au)); 136 CHKERRQ(CheckError(adm, au, user.funcs)); 137 CHKERRQ(MatDestroy(&Interp)); 138 CHKERRQ(DMRestoreGlobalVector(dm, &u)); 139 CHKERRQ(DMDestroy(&dm)); 140 dm = adm; 141 u = au; 142 } 143 CHKERRQ(DMRestoreGlobalVector(dm, &u)); 144 CHKERRQ(DMDestroy(&dm)); 145 ierr = PetscFinalize(); 146 return ierr; 147 } 148 149 /*TEST 150 151 test: 152 suffix: 0 153 args: -num_refine 4 -petscspace_degree 3 \ 154 -dm_plex_dim 1 -dm_plex_box_faces 5 -dm_plex_transform_type refine_1d -dm_plex_hash_location -dm_view 155 156 TEST*/ 157