1 static char help[] = "Tests mesh adaptation with DMPlex and pragmatic.\n"; 2 3 #include <petsc/private/dmpleximpl.h> 4 5 #include <petscksp.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, sizeof(options->mshNam), 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, sizeof(options->bdLabel), 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(); 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 DM cdm, mdm; 82 PetscSection csec, msec; 83 Vec coordinates; 84 const PetscScalar *coords; 85 PetscScalar *met; 86 PetscReal h, *lambda, lbd, lmax; 87 PetscInt pStart, pEnd, p, d; 88 const PetscInt dim = user->dim, Nd = dim*dim; 89 PetscErrorCode ierr; 90 91 PetscFunctionBeginUser; 92 ierr = PetscCalloc1(PetscMax(3, dim),&lambda);CHKERRQ(ierr); 93 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 94 ierr = DMClone(cdm, &mdm);CHKERRQ(ierr); 95 ierr = DMGetLocalSection(cdm, &csec);CHKERRQ(ierr); 96 97 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &msec);CHKERRQ(ierr); 98 ierr = PetscSectionSetNumFields(msec, 1);CHKERRQ(ierr); 99 ierr = PetscSectionSetFieldComponents(msec, 0, Nd);CHKERRQ(ierr); 100 ierr = PetscSectionGetChart(csec, &pStart, &pEnd);CHKERRQ(ierr); 101 ierr = PetscSectionSetChart(msec, pStart, pEnd);CHKERRQ(ierr); 102 for (p = pStart; p < pEnd; ++p) { 103 ierr = PetscSectionSetDof(msec, p, Nd);CHKERRQ(ierr); 104 ierr = PetscSectionSetFieldDof(msec, p, 0, Nd);CHKERRQ(ierr); 105 } 106 ierr = PetscSectionSetUp(msec);CHKERRQ(ierr); 107 ierr = DMSetLocalSection(mdm, msec);CHKERRQ(ierr); 108 ierr = PetscSectionDestroy(&msec);CHKERRQ(ierr); 109 110 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 111 ierr = DMCreateLocalVector(mdm, metric);CHKERRQ(ierr); 112 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 113 ierr = VecGetArray(*metric, &met);CHKERRQ(ierr); 114 for (p = pStart; p < pEnd; ++p) { 115 PetscScalar *pcoords; 116 PetscScalar *pmet; 117 118 ierr = DMPlexPointLocalRead(cdm, p, coords, &pcoords);CHKERRQ(ierr); 119 switch (user->metOpt) { 120 case 0: 121 lbd = 1/(user->hmax*user->hmax); 122 lambda[0] = lambda[1] = lambda[2] = lbd; 123 break; 124 case 1: 125 h = user->hmax - (user->hmax-user->hmin)*PetscRealPart(pcoords[0]); 126 h = h*h; 127 lmax = 1/(user->hmax*user->hmax); 128 lambda[0] = 1/h; 129 lambda[1] = lmax; 130 lambda[2] = lmax; 131 break; 132 case 2: 133 h = user->hmax*PetscAbsReal(((PetscReal) 1.0)-PetscExpReal(-PetscAbsScalar(pcoords[0]-(PetscReal)0.5))) + user->hmin; 134 lbd = 1/(h*h); 135 lmax = 1/(user->hmax*user->hmax); 136 lambda[0] = lbd; 137 lambda[1] = lmax; 138 lambda[2] = lmax; 139 break; 140 default: 141 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "metOpt = 0, 1 or 2, cannot be %d", user->metOpt); 142 } 143 /* Only set the diagonal */ 144 ierr = DMPlexPointLocalRef(mdm, p, met, &pmet);CHKERRQ(ierr); 145 for (d = 0; d < dim; ++d) pmet[d*(dim+1)] = lambda[d]; 146 } 147 ierr = VecRestoreArray(*metric, &met);CHKERRQ(ierr); 148 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 149 ierr = DMDestroy(&mdm);CHKERRQ(ierr); 150 ierr = PetscFree(lambda);CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 154 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 155 { 156 u[0] = x[0] + x[1]; 157 return 0; 158 } 159 160 static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, 161 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 162 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 163 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 164 { 165 g0[0] = 1.0; 166 } 167 168 static PetscErrorCode TestL2Projection(DM dm, DM dma, AppCtx *user) 169 { 170 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *); 171 KSP ksp; 172 PetscDS prob; 173 PetscFE fe; 174 Mat Interp, mass; 175 Vec u, ua, scaling, ones, massLumped, rhs, uproj; 176 PetscReal error; 177 PetscInt dim; 178 MPI_Comm comm; 179 PetscErrorCode ierr; 180 181 PetscFunctionBeginUser; 182 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 183 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 184 ierr = PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, NULL, -1, &fe);CHKERRQ(ierr); 185 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 186 ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe);CHKERRQ(ierr); 187 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 188 ierr = PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, NULL, -1, &fe);CHKERRQ(ierr); 189 ierr = DMGetDS(dma, &prob);CHKERRQ(ierr); 190 ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe);CHKERRQ(ierr); 191 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 192 ierr = PetscDSSetJacobian(prob, 0, 0, identity, NULL, NULL, NULL);CHKERRQ(ierr); 193 194 funcs[0] = linear; 195 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 196 ierr = DMGetGlobalVector(dma, &ua);CHKERRQ(ierr); 197 ierr = DMGetGlobalVector(dma, &ones);CHKERRQ(ierr); 198 ierr = DMGetGlobalVector(dma, &massLumped);CHKERRQ(ierr); 199 ierr = DMGetGlobalVector(dma, &rhs);CHKERRQ(ierr); 200 ierr = DMGetGlobalVector(dma, &uproj);CHKERRQ(ierr); 201 ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, u);CHKERRQ(ierr); 202 ierr = PetscObjectSetName((PetscObject) u, "Original");CHKERRQ(ierr); 203 ierr = VecViewFromOptions(u, NULL, "-orig_vec_view");CHKERRQ(ierr); 204 ierr = DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &error);CHKERRQ(ierr); 205 ierr = PetscPrintf(PETSC_COMM_WORLD, "Original L2 Error: %g\n", (double) error);CHKERRQ(ierr); 206 ierr = DMCreateInterpolation(dm, dma, &Interp, &scaling);CHKERRQ(ierr); 207 ierr = MatInterpolate(Interp, u, ua);CHKERRQ(ierr); 208 ierr = MatDestroy(&Interp);CHKERRQ(ierr); 209 ierr = VecDestroy(&scaling);CHKERRQ(ierr); 210 ierr = PetscObjectSetName((PetscObject) ua, "Interpolation");CHKERRQ(ierr); 211 ierr = VecViewFromOptions(ua, NULL, "-interp_vec_view");CHKERRQ(ierr); 212 ierr = DMComputeL2Diff(dma, 0.0, funcs, NULL, ua, &error);CHKERRQ(ierr); 213 ierr = PetscPrintf(PETSC_COMM_WORLD, "Interpolated L2 Error: %g\n", (double) error);CHKERRQ(ierr); 214 215 ierr = VecSet(ones, 1.0);CHKERRQ(ierr); 216 ierr = DMPlexComputeJacobianAction(dma, NULL, 0, 0, ua, NULL, ones, massLumped, user);CHKERRQ(ierr); 217 ierr = PetscObjectSetName((PetscObject) massLumped, "Lumped mass");CHKERRQ(ierr); 218 ierr = VecViewFromOptions(massLumped, NULL, "-mass_vec_view");CHKERRQ(ierr); 219 ierr = DMCreateMassMatrix(dm, dma, &mass);CHKERRQ(ierr); 220 ierr = MatMult(mass, u, rhs);CHKERRQ(ierr); 221 ierr = MatDestroy(&mass);CHKERRQ(ierr); 222 ierr = VecViewFromOptions(rhs, NULL, "-lumped_rhs_view");CHKERRQ(ierr); 223 ierr = VecPointwiseDivide(uproj, rhs, massLumped);CHKERRQ(ierr); 224 ierr = PetscObjectSetName((PetscObject) uproj, "Different Lumped Projection");CHKERRQ(ierr); 225 ierr = VecViewFromOptions(uproj, NULL, "-lumped_rhs_vec_view");CHKERRQ(ierr); 226 ierr = DMComputeL2Diff(dma, 0.0, funcs, NULL, uproj, &error);CHKERRQ(ierr); 227 ierr = PetscPrintf(PETSC_COMM_WORLD, "Lumped (rhs) L2 Error: %g\n", (double) error);CHKERRQ(ierr); 228 229 ierr = DMCreateMatrix(dma, &mass);CHKERRQ(ierr); 230 ierr = DMPlexSNESComputeJacobianFEM(dma, ua, mass, mass, user);CHKERRQ(ierr); 231 ierr = MatViewFromOptions(mass, NULL, "-mass_mat_view");CHKERRQ(ierr); 232 ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);CHKERRQ(ierr); 233 ierr = KSPSetOperators(ksp, mass, mass);CHKERRQ(ierr); 234 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 235 ierr = KSPSolve(ksp, rhs, uproj);CHKERRQ(ierr); 236 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 237 ierr = PetscObjectSetName((PetscObject) uproj, "Full Projection");CHKERRQ(ierr); 238 ierr = VecViewFromOptions(uproj, NULL, "-proj_vec_view");CHKERRQ(ierr); 239 ierr = DMComputeL2Diff(dma, 0.0, funcs, NULL, uproj, &error);CHKERRQ(ierr); 240 ierr = PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double) error);CHKERRQ(ierr); 241 ierr = MatDestroy(&mass);CHKERRQ(ierr); 242 243 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 244 ierr = DMRestoreGlobalVector(dma, &ua);CHKERRQ(ierr); 245 ierr = DMRestoreGlobalVector(dma, &ones);CHKERRQ(ierr); 246 ierr = DMRestoreGlobalVector(dma, &massLumped);CHKERRQ(ierr); 247 ierr = DMRestoreGlobalVector(dma, &rhs);CHKERRQ(ierr); 248 ierr = DMRestoreGlobalVector(dma, &uproj);CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 252 int main (int argc, char * argv[]) { 253 AppCtx user; /* user-defined work context */ 254 DMLabel bdLabel = NULL; 255 MPI_Comm comm; 256 DM dma, odm; 257 Vec metric; 258 size_t len; 259 PetscErrorCode ierr; 260 261 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 262 comm = PETSC_COMM_WORLD; 263 ierr = ProcessOptions(comm, &user);CHKERRQ(ierr); 264 265 ierr = CreateMesh(comm, &user);CHKERRQ(ierr); 266 ierr = PetscObjectSetName((PetscObject) user.dm, "DMinit");CHKERRQ(ierr); 267 ierr = DMViewFromOptions(user.dm, NULL, "-init_dm_view");CHKERRQ(ierr); 268 269 odm = user.dm; 270 ierr = DMPlexDistributeOverlap(odm, 1, NULL, &user.dm);CHKERRQ(ierr); 271 if (!user.dm) {user.dm = odm;} 272 else {ierr = DMDestroy(&odm);CHKERRQ(ierr);} 273 ierr = ComputeMetric(user.dm, &user, &metric);CHKERRQ(ierr); 274 ierr = PetscStrlen(user.bdLabel, &len);CHKERRQ(ierr); 275 if (len) { 276 ierr = DMCreateLabel(user.dm, user.bdLabel);CHKERRQ(ierr); 277 ierr = DMGetLabel(user.dm, user.bdLabel, &bdLabel);CHKERRQ(ierr); 278 } 279 ierr = DMAdaptMetric(user.dm, metric, bdLabel, &dma);CHKERRQ(ierr); 280 ierr = PetscObjectSetName((PetscObject) dma, "DMadapt");CHKERRQ(ierr); 281 ierr = PetscObjectSetOptionsPrefix((PetscObject) dma, "adapt_");CHKERRQ(ierr); 282 ierr = DMViewFromOptions(dma, NULL, "-dm_view");CHKERRQ(ierr); 283 if (user.doL2) {ierr = TestL2Projection(user.dm, dma, &user);CHKERRQ(ierr);} 284 ierr = DMDestroy(&dma);CHKERRQ(ierr); 285 ierr = VecDestroy(&metric);CHKERRQ(ierr); 286 ierr = DMDestroy(&user.dm);CHKERRQ(ierr); 287 PetscFinalize(); 288 return 0; 289 } 290 291 /*TEST 292 293 test: 294 suffix: 0 295 requires: pragmatic 296 TODO: broken 297 args: -dim 2 -nbrVerEdge 5 -dm_plex_separate_marker 0 -met 2 -init_dm_view -adapt_dm_view 298 test: 299 suffix: 1 300 requires: pragmatic 301 TODO: broken 302 args: -dim 2 -nbrVerEdge 5 -dm_plex_separate_marker 1 -bdLabel marker -met 2 -init_dm_view -adapt_dm_view 303 test: 304 suffix: 2 305 requires: pragmatic 306 TODO: broken 307 args: -dim 3 -nbrVerEdge 5 -met 2 -init_dm_view -adapt_dm_view 308 test: 309 suffix: 3 310 requires: pragmatic 311 TODO: broken 312 args: -dim 3 -nbrVerEdge 5 -bdLabel marker -met 2 -init_dm_view -adapt_dm_view 313 test: 314 suffix: 4 315 requires: pragmatic 316 TODO: broken 317 nsize: 2 318 args: -dim 2 -nbrVerEdge 3 -dm_plex_separate_marker 0 -met 2 -init_dm_view -adapt_dm_view 319 test: 320 suffix: 5 321 requires: pragmatic 322 TODO: broken 323 nsize: 4 324 args: -dim 2 -nbrVerEdge 3 -dm_plex_separate_marker 0 -met 2 -init_dm_view -adapt_dm_view 325 test: 326 suffix: 6 327 requires: pragmatic 328 TODO: broken 329 nsize: 2 330 args: -dim 3 -nbrVerEdge 10 -dm_plex_separate_marker 0 -met 0 -hmin 0.01 -hmax 0.03 -init_dm_view -adapt_dm_view 331 test: 332 suffix: 7 333 requires: pragmatic 334 TODO: broken 335 nsize: 5 336 args: -dim 2 -nbrVerEdge 20 -dm_plex_separate_marker 0 -met 2 -hmax 0.5 -hmin 0.001 -init_dm_view -adapt_dm_view 337 test: 338 suffix: proj_0 339 requires: pragmatic 340 TODO: broken 341 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 342 test: 343 suffix: proj_1 344 requires: pragmatic 345 TODO: broken 346 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 347 348 TEST*/ 349