1 static char help[] = "Tests dof numberings for external integrators such as LibCEED.\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petscds.h> 5 6 typedef struct { 7 PetscBool useFE; 8 PetscInt check_face; 9 PetscBool closure_tensor; 10 } AppCtx; 11 12 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 13 { 14 PetscFunctionBeginUser; 15 options->useFE = PETSC_TRUE; 16 options->check_face = 1; 17 options->closure_tensor = PETSC_FALSE; 18 PetscOptionsBegin(comm, "", "Dof Ordering Options", "DMPLEX"); 19 PetscCall(PetscOptionsBool("-use_fe", "Use FE or FV discretization", "ex49.c", options->useFE, &options->useFE, NULL)); 20 PetscCall(PetscOptionsInt("-check_face", "Face set to report on", "ex49.c", options->check_face, &options->check_face, NULL)); 21 PetscCall(PetscOptionsBool("-closure_tensor", "Use DMPlexSetClosurePermutationTensor()", "ex49.c", options->closure_tensor, &options->closure_tensor, NULL)); 22 PetscOptionsEnd(); 23 PetscFunctionReturn(PETSC_SUCCESS); 24 } 25 26 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 27 { 28 PetscFunctionBeginUser; 29 PetscCall(DMCreate(comm, dm)); 30 PetscCall(DMSetType(*dm, DMPLEX)); 31 PetscCall(DMSetFromOptions(*dm)); 32 PetscCall(DMSetApplicationContext(*dm, user)); 33 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 34 PetscFunctionReturn(PETSC_SUCCESS); 35 } 36 37 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 38 { 39 DM cdm = dm; 40 PetscInt dim; 41 42 PetscFunctionBeginUser; 43 PetscCall(DMGetDimension(dm, &dim)); 44 if (user->useFE) { 45 PetscFE fe; 46 DMPolytopeType ct; 47 PetscInt cStart; 48 49 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 50 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 51 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, NULL, -1, &fe)); 52 PetscCall(PetscObjectSetName((PetscObject)fe, "scalar")); 53 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 54 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe)); 55 PetscCall(PetscFEDestroy(&fe)); 56 } else { 57 PetscFV fv; 58 59 PetscCall(PetscFVCreate(PETSC_COMM_SELF, &fv)); 60 PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES)); 61 PetscCall(PetscFVSetNumComponents(fv, dim)); 62 PetscCall(PetscFVSetSpatialDimension(fv, dim)); 63 PetscCall(PetscFVSetFromOptions(fv)); 64 PetscCall(PetscFVSetUp(fv)); 65 PetscCall(PetscObjectSetName((PetscObject)fv, "vector")); 66 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fv)); 67 PetscCall(PetscFVDestroy(&fv)); 68 } 69 PetscCall(DMCreateDS(dm)); 70 while (cdm) { 71 PetscCall(DMCopyDisc(dm, cdm)); 72 PetscCall(DMGetCoarseDM(cdm, &cdm)); 73 } 74 PetscFunctionReturn(PETSC_SUCCESS); 75 } 76 77 static PetscErrorCode CheckOffsets(DM dm, AppCtx *user, const char *domain_name, PetscInt label_value, PetscInt height) 78 { 79 const char *height_name[] = {"cells", "faces"}; 80 DMLabel domain_label = NULL; 81 DM cdm; 82 PetscInt Nf, f; 83 ISLocalToGlobalMapping ltog; 84 85 PetscFunctionBeginUser; 86 if (domain_name) PetscCall(DMGetLabel(dm, domain_name, &domain_label)); 87 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "## %s: '%s' {%" PetscInt_FMT "}%s\n", height_name[height], domain_name ? domain_name : "default", label_value, domain_name && !domain_label ? " (null label)" : "")); 88 if (domain_name && !domain_label) PetscFunctionReturn(PETSC_SUCCESS); 89 if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 90 // Offsets for cell closures 91 PetscCall(DMGetNumFields(dm, &Nf)); 92 for (f = 0; f < Nf; ++f) { 93 PetscObject obj; 94 PetscClassId id; 95 char name[PETSC_MAX_PATH_LEN]; 96 97 PetscCall(DMGetField(dm, f, NULL, &obj)); 98 PetscCall(PetscObjectGetClassId(obj, &id)); 99 if (id == PETSCFE_CLASSID) { 100 IS offIS; 101 PetscInt *offsets, Ncell, Ncl, Nc, n; 102 103 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets)); 104 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS)); 105 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f)); 106 PetscCall(PetscObjectSetName((PetscObject)offIS, name)); 107 PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view")); 108 PetscCall(ISDestroy(&offIS)); 109 } else if (id == PETSCFV_CLASSID) { 110 IS offIS; 111 PetscInt *offsets, *offsetsNeg, *offsetsPos, Nface, Nc, n, i = 0; 112 113 PetscCall(DMPlexGetLocalOffsetsSupport(dm, domain_label, label_value, &Nface, &Nc, &n, &offsetsNeg, &offsetsPos)); 114 PetscCall(PetscMalloc1(Nface * Nc * 2, &offsets)); 115 for (PetscInt f = 0; f < Nface; ++f) { 116 for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsNeg[f] + c; 117 for (PetscInt c = 0; c < Nc; ++c) offsets[i++] = offsetsPos[f] + c; 118 } 119 PetscCheck(i == Nface * Nc * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total offsets %" PetscInt_FMT " != %" PetscInt_FMT, i, Nface * Nc * 2); 120 PetscCall(PetscFree(offsetsNeg)); 121 PetscCall(PetscFree(offsetsPos)); 122 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Nface * Nc * 2, offsets, PETSC_OWN_POINTER, &offIS)); 123 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f)); 124 PetscCall(PetscObjectSetName((PetscObject)offIS, name)); 125 PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view")); 126 PetscCall(ISDestroy(&offIS)); 127 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unrecognized type for DM field %" PetscInt_FMT, f); 128 } 129 PetscCall(DMGetLocalToGlobalMapping(dm, <og)); 130 PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-ltog_view")); 131 132 // Offsets for coordinates 133 { 134 Vec X; 135 PetscSection s; 136 const PetscScalar *x; 137 const char *cname; 138 PetscInt cdim, *offsets, Ncell, Ncl, Nc, n; 139 PetscBool isDG = PETSC_FALSE; 140 141 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 142 if (!cdm) { 143 PetscCall(DMGetCoordinateDM(dm, &cdm)); 144 cname = "Coordinates"; 145 PetscCall(DMGetCoordinatesLocal(dm, &X)); 146 } else { 147 isDG = PETSC_TRUE; 148 cname = "DG Coordinates"; 149 PetscCall(DMGetCellCoordinatesLocal(dm, &X)); 150 } 151 if (isDG && height) PetscFunctionReturn(PETSC_SUCCESS); 152 if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label)); 153 if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL)); 154 PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets)); 155 PetscCall(DMGetCoordinateDim(dm, &cdim)); 156 PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim); 157 PetscCall(DMGetLocalSection(cdm, &s)); 158 PetscCall(VecGetArrayRead(X, &x)); 159 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element in %s order\n", cname, user->closure_tensor ? "tensor" : "bfs")); 160 for (PetscInt c = 0; c < Ncell; ++c) { 161 for (PetscInt v = 0; v < Ncl; ++v) { 162 PetscInt off = offsets[c * Ncl + v], dgdof; 163 const PetscScalar *vx = &x[off]; 164 165 if (isDG) { 166 PetscCall(PetscSectionGetDof(s, c, &dgdof)); 167 PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof); 168 } 169 switch (cdim) { 170 case 1: 171 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]))); 172 break; 173 case 2: 174 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1]))); 175 break; 176 case 3: 177 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f, % 4.2f, % 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]), (double)PetscRealPart(vx[1]), (double)PetscRealPart(vx[2]))); 178 } 179 } 180 } 181 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout)); 182 PetscCall(VecRestoreArrayRead(X, &x)); 183 PetscCall(PetscFree(offsets)); 184 PetscCall(DMGetLocalToGlobalMapping(cdm, <og)); 185 PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-coord_ltog_view")); 186 { 187 DM clonedm; 188 Vec cloneX, X; 189 PetscInt clone_num_x, num_x; 190 const PetscScalar *clonex, *x; 191 192 PetscCall(DMClone(dm, &clonedm)); 193 { // Force recreation of local coordinate vector 194 Vec X_global; 195 196 PetscCall(DMGetCoordinates(dm, &X_global)); 197 PetscCall(DMSetCoordinates(clonedm, X_global)); 198 } 199 PetscCall(DMGetCoordinatesLocal(dm, &X)); 200 PetscCall(DMGetCoordinatesLocal(clonedm, &cloneX)); 201 PetscCall(VecGetLocalSize(X, &num_x)); 202 PetscCall(VecGetLocalSize(cloneX, &clone_num_x)); 203 PetscCheck(num_x == clone_num_x, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Cloned DM coordinate size (%" PetscInt_FMT ") different from original DM coordinate size (%" PetscInt_FMT ")", clone_num_x, num_x); 204 205 PetscCall(VecGetArrayRead(X, &x)); 206 PetscCall(VecGetArrayRead(cloneX, &clonex)); 207 208 for (PetscInt i = 0; i < num_x; i++) { 209 PetscCheck(PetscIsCloseAtTolScalar(x[i], clonex[i], 1e-13, 1e-13), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Original coordinate (%4.2f) and cloned coordinate (%4.2f) are different", (double)PetscRealPart(x[i]), (double)PetscRealPart(clonex[i])); 210 } 211 212 PetscCall(VecRestoreArrayRead(X, &x)); 213 PetscCall(VecRestoreArrayRead(cloneX, &clonex)); 214 PetscCall(DMDestroy(&clonedm)); 215 } 216 } 217 PetscFunctionReturn(PETSC_SUCCESS); 218 } 219 220 int main(int argc, char **argv) 221 { 222 DM dm; 223 AppCtx user; 224 PetscInt depth; 225 226 PetscFunctionBeginUser; 227 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 228 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 229 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 230 PetscCall(SetupDiscretization(dm, &user)); 231 PetscCall(CheckOffsets(dm, &user, NULL, 0, 0)); 232 PetscCall(DMPlexGetDepth(dm, &depth)); 233 if (depth > 1) PetscCall(CheckOffsets(dm, &user, "Face Sets", user.check_face, 1)); 234 PetscCall(DMDestroy(&dm)); 235 PetscCall(PetscFinalize()); 236 return 0; 237 } 238 239 /*TEST 240 241 test: 242 suffix: 0 243 requires: triangle 244 args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view 245 246 test: 247 suffix: 1 248 args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \ 249 -dm_view -offsets_view 250 251 test: 252 suffix: cg_2d 253 args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \ 254 -dm_view -offsets_view 255 256 test: 257 suffix: 1d_sfc 258 args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view 259 260 test: 261 suffix: 2d_sfc 262 nsize: 2 263 args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 4,3 -dm_distribute 0 -petscspace_degree 1 -dm_view 264 265 test: 266 suffix: 2d_sfc_periodic 267 nsize: 2 268 args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 4,3 -dm_distribute 0 -petscspace_degree 1 -dm_plex_box_bd periodic,none -dm_view ::ascii_info_detail 269 270 testset: 271 args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 3,2 -petscspace_degree 1 -dm_view ::ascii_info_detail -closure_tensor 272 nsize: 2 273 test: 274 suffix: 2d_sfc_periodic_stranded 275 args: -dm_distribute 0 -dm_plex_box_bd none,periodic 276 test: 277 suffix: 2d_sfc_periodic_stranded_dist 278 args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd none,periodic 279 test: 280 suffix: 2d_sfc_biperiodic_stranded 281 args: -dm_distribute 0 -dm_plex_box_bd periodic,periodic 282 test: 283 suffix: 2d_sfc_biperiodic_stranded_dist 284 args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_bd periodic,periodic 285 test: 286 suffix: 2d_sfc_biperiodic_stranded_dist_box_label 287 args: -dm_distribute 1 -petscpartitioner_type simple -dm_plex_box_label_bd periodic,periodic -dm_plex_box_label 288 289 test: 290 suffix: fv_0 291 requires: triangle 292 args: -dm_refine 1 -use_fe 0 -dm_view -offsets_view 293 294 TEST*/ 295