1 static char help[] = "Tests mesh adaptation with DMPlex and pragmatic.\n"; 2 3 #include <petsc/private/dmpleximpl.h> 4 5 #include <petscsnes.h> 6 7 typedef struct { 8 DM dm; 9 /* Definition of the test case (mesh and metric field) */ 10 PetscInt dim; /* The topological mesh dimension */ 11 char mshNam[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */ 12 PetscInt nbrVerEdge; /* Number of vertices per edge if unit square/cube generated */ 13 char bdLabel[PETSC_MAX_PATH_LEN]; /* Name of the label marking boundary facets */ 14 PetscInt metOpt; /* Different choices of metric */ 15 PetscReal hmax, hmin; /* Max and min sizes prescribed by the metric */ 16 PetscBool doL2; /* Test L2 projection */ 17 } AppCtx; 18 19 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 20 { 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 options->dim = 2; 25 ierr = PetscStrcpy(options->mshNam, "");CHKERRQ(ierr); 26 options->nbrVerEdge = 5; 27 ierr = PetscStrcpy(options->bdLabel, "");CHKERRQ(ierr); 28 options->metOpt = 1; 29 options->hmin = 0.05; 30 options->hmax = 0.5; 31 options->doL2 = PETSC_FALSE; 32 33 ierr = PetscOptionsBegin(comm, "", "Meshing Adaptation Options", "DMPLEX");CHKERRQ(ierr); 34 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex19.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 35 ierr = PetscOptionsString("-msh", "Name of the mesh filename if any", "ex19.c", options->mshNam, options->mshNam, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 36 ierr = PetscOptionsBoundedInt("-nbrVerEdge", "Number of vertices per edge if unit square/cube generated", "ex19.c", options->nbrVerEdge, &options->nbrVerEdge, NULL,0);CHKERRQ(ierr); 37 ierr = PetscOptionsString("-bdLabel", "Name of the label marking boundary facets", "ex19.c", options->bdLabel, options->bdLabel, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 38 ierr = PetscOptionsBoundedInt("-met", "Different choices of metric", "ex19.c", options->metOpt, &options->metOpt, NULL,0);CHKERRQ(ierr); 39 ierr = PetscOptionsReal("-hmax", "Max size prescribed by the metric", "ex19.c", options->hmax, &options->hmax, NULL);CHKERRQ(ierr); 40 ierr = PetscOptionsReal("-hmin", "Min size prescribed by the metric", "ex19.c", options->hmin, &options->hmin, NULL);CHKERRQ(ierr); 41 ierr = PetscOptionsBool("-do_L2", "Test L2 projection", "ex19.c", options->doL2, &options->doL2, NULL);CHKERRQ(ierr); 42 ierr = PetscOptionsEnd();CHKERRQ(ierr); 43 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user) 48 { 49 PetscBool flag; 50 PetscErrorCode ierr; 51 52 PetscFunctionBegin; 53 ierr = PetscStrcmp(user->mshNam, "", &flag);CHKERRQ(ierr); 54 if (flag) { 55 PetscInt faces[3]; 56 faces[0] = user->nbrVerEdge-1; 57 faces[1] = user->nbrVerEdge-1; 58 faces[2] = user->nbrVerEdge-1; 59 60 ierr = DMPlexCreateBoxMesh(comm, user->dim, PETSC_TRUE, faces, NULL, NULL, NULL, PETSC_TRUE, &user->dm);CHKERRQ(ierr); 61 } else { 62 ierr = DMPlexCreateFromFile(comm, user->mshNam, PETSC_TRUE, &user->dm);CHKERRQ(ierr); 63 ierr = DMGetDimension(user->dm, &user->dim);CHKERRQ(ierr); 64 } 65 { 66 DM distributedMesh = NULL; 67 68 /* Distribute mesh over processes */ 69 ierr = DMPlexDistribute(user->dm, 0, NULL, &distributedMesh);CHKERRQ(ierr); 70 if (distributedMesh) { 71 ierr = DMDestroy(&user->dm);CHKERRQ(ierr); 72 user->dm = distributedMesh; 73 } 74 } 75 ierr = DMSetFromOptions(user->dm);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode ComputeMetric(DM dm, AppCtx *user, Vec *metric) 80 { 81 PetscReal lambda; 82 PetscErrorCode ierr; 83 84 PetscFunctionBeginUser; 85 86 /* Specify a uniform, isotropic metric */ 87 lambda = 1/(user->hmax*user->hmax); 88 ierr = DMPlexMetricCreateUniform(dm, 0, lambda, metric);CHKERRQ(ierr); 89 90 /* Modify the first diagonal entry, if requested */ 91 if (user->metOpt != 0) { 92 DM cdm; 93 PetscSection csec; 94 Vec coordinates; 95 const PetscScalar *coords; 96 PetscScalar *met; 97 PetscReal h; 98 PetscInt pStart, pEnd, p; 99 100 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 101 ierr = DMGetLocalSection(cdm, &csec);CHKERRQ(ierr); 102 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 103 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 104 ierr = VecGetArray(*metric, &met);CHKERRQ(ierr); 105 ierr = PetscSectionGetChart(csec, &pStart, &pEnd);CHKERRQ(ierr); 106 for (p = pStart; p < pEnd; ++p) { 107 PetscScalar *pcoords; 108 PetscScalar *pmet; 109 110 ierr = DMPlexPointLocalRead(cdm, p, coords, &pcoords);CHKERRQ(ierr); 111 switch (user->metOpt) { 112 case 1: 113 h = user->hmax - (user->hmax-user->hmin)*PetscRealPart(pcoords[0]); 114 break; 115 case 2: 116 h = user->hmax*PetscAbsReal(((PetscReal) 1.0)-PetscExpReal(-PetscAbsScalar(pcoords[0]-(PetscReal)0.5))) + user->hmin; 117 break; 118 default: 119 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "metOpt = 0, 1 or 2, cannot be %d", user->metOpt); 120 } 121 lambda = 1/(h*h); 122 ierr = DMPlexPointLocalRef(dm, p, met, &pmet);CHKERRQ(ierr); 123 pmet[0] = lambda; 124 } 125 ierr = VecRestoreArray(*metric, &met);CHKERRQ(ierr); 126 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 127 } 128 PetscFunctionReturn(0); 129 } 130 131 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 132 { 133 u[0] = x[0] + x[1]; 134 return 0; 135 } 136 137 static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, 138 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 139 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 140 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 141 { 142 g0[0] = 1.0; 143 } 144 145 static PetscErrorCode TestL2Projection(DM dm, DM dma, AppCtx *user) 146 { 147 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *); 148 KSP ksp; 149 PetscDS prob; 150 PetscFE fe; 151 Mat Interp, mass; 152 Vec u, ua, scaling, ones, massLumped, rhs, uproj; 153 PetscReal error; 154 PetscInt dim; 155 MPI_Comm comm; 156 PetscErrorCode ierr; 157 158 PetscFunctionBeginUser; 159 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 160 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 161 ierr = PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, NULL, -1, &fe);CHKERRQ(ierr); 162 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 163 ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe);CHKERRQ(ierr); 164 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 165 ierr = PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, NULL, -1, &fe);CHKERRQ(ierr); 166 ierr = DMGetDS(dma, &prob);CHKERRQ(ierr); 167 ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe);CHKERRQ(ierr); 168 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 169 ierr = PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL);CHKERRQ(ierr); 170 171 funcs[0] = linear; 172 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 173 ierr = DMGetGlobalVector(dma, &ua);CHKERRQ(ierr); 174 ierr = DMGetGlobalVector(dma, &ones);CHKERRQ(ierr); 175 ierr = DMGetGlobalVector(dma, &massLumped);CHKERRQ(ierr); 176 ierr = DMGetGlobalVector(dma, &rhs);CHKERRQ(ierr); 177 ierr = DMGetGlobalVector(dma, &uproj);CHKERRQ(ierr); 178 ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, u);CHKERRQ(ierr); 179 ierr = PetscObjectSetName((PetscObject) u, "Original");CHKERRQ(ierr); 180 ierr = VecViewFromOptions(u, NULL, "-orig_vec_view");CHKERRQ(ierr); 181 ierr = DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error);CHKERRQ(ierr); 182 ierr = PetscPrintf(PETSC_COMM_WORLD, "Original L2 Error: %g\n", (double) error);CHKERRQ(ierr); 183 ierr = DMCreateInterpolation(dm, dma, &Interp, &scaling);CHKERRQ(ierr); 184 ierr = MatInterpolate(Interp, u, ua);CHKERRQ(ierr); 185 ierr = MatDestroy(&Interp);CHKERRQ(ierr); 186 ierr = VecDestroy(&scaling);CHKERRQ(ierr); 187 ierr = PetscObjectSetName((PetscObject) ua, "Interpolation");CHKERRQ(ierr); 188 ierr = VecViewFromOptions(ua, NULL, "-interp_vec_view");CHKERRQ(ierr); 189 ierr = DMComputeL2Diff(dma, 0.0, funcs, NULL, ua, &error);CHKERRQ(ierr); 190 ierr = PetscPrintf(PETSC_COMM_WORLD, "Interpolated L2 Error: %g\n", (double) error);CHKERRQ(ierr); 191 192 ierr = VecSet(ones, 1.0);CHKERRQ(ierr); 193 ierr = DMSNESComputeJacobianAction(dma, ua, ones, massLumped, user);CHKERRQ(ierr); 194 ierr = PetscObjectSetName((PetscObject) massLumped, "Lumped mass");CHKERRQ(ierr); 195 ierr = VecViewFromOptions(massLumped, NULL, "-mass_vec_view");CHKERRQ(ierr); 196 ierr = DMCreateMassMatrix(dm, dma, &mass);CHKERRQ(ierr); 197 ierr = MatMult(mass, u, rhs);CHKERRQ(ierr); 198 ierr = MatDestroy(&mass);CHKERRQ(ierr); 199 ierr = VecViewFromOptions(rhs, NULL, "-lumped_rhs_view");CHKERRQ(ierr); 200 ierr = VecPointwiseDivide(uproj, rhs, massLumped);CHKERRQ(ierr); 201 ierr = PetscObjectSetName((PetscObject) uproj, "Different Lumped Projection");CHKERRQ(ierr); 202 ierr = VecViewFromOptions(uproj, NULL, "-lumped_rhs_vec_view");CHKERRQ(ierr); 203 ierr = DMComputeL2Diff(dma, 0.0, funcs, NULL, uproj, &error);CHKERRQ(ierr); 204 ierr = PetscPrintf(PETSC_COMM_WORLD, "Lumped (rhs) L2 Error: %g\n", (double) error);CHKERRQ(ierr); 205 206 ierr = DMCreateMatrix(dma, &mass);CHKERRQ(ierr); 207 ierr = DMPlexSNESComputeJacobianFEM(dma, ua, mass, mass, user);CHKERRQ(ierr); 208 ierr = MatViewFromOptions(mass, NULL, "-mass_mat_view");CHKERRQ(ierr); 209 ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);CHKERRQ(ierr); 210 ierr = KSPSetOperators(ksp, mass, mass);CHKERRQ(ierr); 211 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 212 ierr = KSPSolve(ksp, rhs, uproj);CHKERRQ(ierr); 213 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 214 ierr = PetscObjectSetName((PetscObject) uproj, "Full Projection");CHKERRQ(ierr); 215 ierr = VecViewFromOptions(uproj, NULL, "-proj_vec_view");CHKERRQ(ierr); 216 ierr = DMComputeL2Diff(dma, 0.0, funcs, NULL, uproj, &error);CHKERRQ(ierr); 217 ierr = PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double) error);CHKERRQ(ierr); 218 ierr = MatDestroy(&mass);CHKERRQ(ierr); 219 220 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 221 ierr = DMRestoreGlobalVector(dma, &ua);CHKERRQ(ierr); 222 ierr = DMRestoreGlobalVector(dma, &ones);CHKERRQ(ierr); 223 ierr = DMRestoreGlobalVector(dma, &massLumped);CHKERRQ(ierr); 224 ierr = DMRestoreGlobalVector(dma, &rhs);CHKERRQ(ierr); 225 ierr = DMRestoreGlobalVector(dma, &uproj);CHKERRQ(ierr); 226 PetscFunctionReturn(0); 227 } 228 229 int main (int argc, char * argv[]) { 230 AppCtx user; /* user-defined work context */ 231 DMLabel bdLabel = NULL; 232 MPI_Comm comm; 233 DM dma, odm; 234 Vec metric; 235 size_t len; 236 PetscErrorCode ierr; 237 238 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 239 comm = PETSC_COMM_WORLD; 240 ierr = ProcessOptions(comm, &user);CHKERRQ(ierr); 241 242 ierr = CreateMesh(comm, &user);CHKERRQ(ierr); 243 ierr = PetscObjectSetName((PetscObject) user.dm, "DMinit");CHKERRQ(ierr); 244 ierr = DMViewFromOptions(user.dm, NULL, "-init_dm_view");CHKERRQ(ierr); 245 246 odm = user.dm; 247 ierr = DMPlexDistributeOverlap(odm, 1, NULL, &user.dm);CHKERRQ(ierr); 248 if (!user.dm) {user.dm = odm;} 249 else {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 250 ierr = ComputeMetric(user.dm, &user, &metric);CHKERRQ(ierr); 251 ierr = PetscStrlen(user.bdLabel, &len);CHKERRQ(ierr); 252 if (len) { 253 ierr = DMCreateLabel(user.dm, user.bdLabel);CHKERRQ(ierr); 254 ierr = DMGetLabel(user.dm, user.bdLabel, &bdLabel);CHKERRQ(ierr); 255 } 256 ierr = DMAdaptMetric(user.dm, metric, bdLabel, &dma);CHKERRQ(ierr); 257 ierr = PetscObjectSetName((PetscObject) dma, "DMadapt");CHKERRQ(ierr); 258 ierr = PetscObjectSetOptionsPrefix((PetscObject) dma, "adapt_");CHKERRQ(ierr); 259 ierr = DMViewFromOptions(dma, NULL, "-dm_view");CHKERRQ(ierr); 260 if (user.doL2) {ierr = TestL2Projection(user.dm, dma, &user);CHKERRQ(ierr);} 261 ierr = DMDestroy(&dma);CHKERRQ(ierr); 262 ierr = VecDestroy(&metric);CHKERRQ(ierr); 263 ierr = DMDestroy(&user.dm);CHKERRQ(ierr); 264 PetscFinalize(); 265 return 0; 266 } 267 268 /*TEST 269 270 build: 271 requires: pragmatic 272 273 test: 274 suffix: 0 275 args: -dim 2 -nbrVerEdge 5 -dm_plex_separate_marker 0 -met 2 -init_dm_view -adapt_dm_view 276 test: 277 suffix: 1 278 args: -dim 2 -nbrVerEdge 5 -dm_plex_separate_marker 1 -bdLabel marker -met 2 -init_dm_view -adapt_dm_view 279 test: 280 suffix: 2 281 args: -dim 3 -nbrVerEdge 5 -met 2 -init_dm_view -adapt_dm_view 282 test: 283 suffix: 3 284 args: -dim 3 -nbrVerEdge 5 -bdLabel marker -met 2 -init_dm_view -adapt_dm_view 285 test: 286 suffix: 4 287 nsize: 2 288 args: -dim 2 -nbrVerEdge 3 -dm_plex_separate_marker 0 -met 2 -init_dm_view -adapt_dm_view 289 test: 290 suffix: 5 291 nsize: 4 292 args: -dim 2 -nbrVerEdge 3 -dm_plex_separate_marker 0 -met 2 -init_dm_view -adapt_dm_view 293 test: 294 suffix: 6 295 nsize: 2 296 args: -dim 3 -nbrVerEdge 10 -met 0 -hmin 0.01 -hmax 0.03 -init_dm_view -adapt_dm_view 297 test: 298 suffix: 7 299 nsize: 2 300 args: -dim 2 -nbrVerEdge 20 -dm_plex_separate_marker 0 -met 2 -hmax 0.5 -hmin 0.001 -init_dm_view -adapt_dm_view 301 test: 302 suffix: proj_0 303 TODO: broken 304 args: -dim 2 -nbrVerEdge 3 -dm_plex_separate_marker 0 -init_dm_view -adapt_dm_view -do_L2 -petscspace_degree 1 -petscfe_default_quadrature_order 1 -dm_plex_hash_location -pc_type lu 305 test: 306 suffix: proj_1 307 TODO: broken 308 args: -dim 2 -nbrVerEdge 5 -dm_plex_separate_marker 0 -init_dm_view -adapt_dm_view -do_L2 -petscspace_degree 2 -petscfe_default_quadrature_order 4 -dm_plex_hash_location -pc_type lu 309 310 TEST*/ 311