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