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