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 PetscInt check_face; 8 PetscBool closure_tensor; 9 } AppCtx; 10 11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 12 { 13 PetscFunctionBeginUser; 14 options->check_face = 1; 15 options->closure_tensor = PETSC_FALSE; 16 PetscOptionsBegin(comm, "", "Dof Ordering Options", "DMPLEX"); 17 PetscCall(PetscOptionsInt("-check_face", "Face set to report on", "ex49.c", options->check_face, &options->check_face, NULL)); 18 PetscCall(PetscOptionsBool("-closure_tensor", "Use DMPlexSetClosurePermutationTensor()", "ex49.c", options->closure_tensor, &options->closure_tensor, NULL)); 19 PetscOptionsEnd(); 20 PetscFunctionReturn(PETSC_SUCCESS); 21 } 22 23 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 24 { 25 PetscFunctionBeginUser; 26 PetscCall(DMCreate(comm, dm)); 27 PetscCall(DMSetType(*dm, DMPLEX)); 28 PetscCall(DMSetFromOptions(*dm)); 29 PetscCall(DMSetApplicationContext(*dm, user)); 30 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 31 PetscFunctionReturn(PETSC_SUCCESS); 32 } 33 34 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 35 { 36 DM cdm = dm; 37 PetscFE fe; 38 PetscInt dim; 39 40 PetscFunctionBeginUser; 41 PetscCall(DMGetDimension(dm, &dim)); 42 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, -1, &fe)); 43 PetscCall(PetscObjectSetName((PetscObject)fe, "scalar")); 44 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 45 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe)); 46 PetscCall(PetscFEDestroy(&fe)); 47 PetscCall(DMCreateDS(dm)); 48 while (cdm) { 49 PetscCall(DMCopyDisc(dm, cdm)); 50 PetscCall(DMGetCoarseDM(cdm, &cdm)); 51 } 52 PetscFunctionReturn(PETSC_SUCCESS); 53 } 54 55 static PetscErrorCode CheckOffsets(DM dm, AppCtx *user, const char *domain_name, PetscInt label_value, PetscInt height) 56 { 57 const char *height_name[] = {"cells", "faces"}; 58 DMLabel domain_label = NULL; 59 DM cdm; 60 IS offIS; 61 PetscInt *offsets, Ncell, Ncl, Nc, n; 62 PetscInt Nf, f; 63 ISLocalToGlobalMapping ltog; 64 65 PetscFunctionBeginUser; 66 if (domain_name) PetscCall(DMGetLabel(dm, domain_name, &domain_label)); 67 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)" : "")); 68 if (domain_name && !domain_label) PetscFunctionReturn(PETSC_SUCCESS); 69 if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 70 // Offsets for cell closures 71 PetscCall(DMGetNumFields(dm, &Nf)); 72 for (f = 0; f < Nf; ++f) { 73 char name[PETSC_MAX_PATH_LEN]; 74 75 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets)); 76 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS)); 77 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f)); 78 PetscCall(PetscObjectSetName((PetscObject)offIS, name)); 79 PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view")); 80 PetscCall(ISDestroy(&offIS)); 81 } 82 PetscCall(DMGetLocalToGlobalMapping(dm, <og)); 83 PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-ltog_view")); 84 85 // Offsets for coordinates 86 { 87 Vec X; 88 PetscSection s; 89 const PetscScalar *x; 90 const char *cname; 91 PetscInt cdim; 92 PetscBool isDG = PETSC_FALSE; 93 94 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 95 if (!cdm) { 96 PetscCall(DMGetCoordinateDM(dm, &cdm)); 97 cname = "Coordinates"; 98 PetscCall(DMGetCoordinatesLocal(dm, &X)); 99 } else { 100 isDG = PETSC_TRUE; 101 cname = "DG Coordinates"; 102 PetscCall(DMGetCellCoordinatesLocal(dm, &X)); 103 } 104 if (isDG && height) PetscFunctionReturn(PETSC_SUCCESS); 105 if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label)); 106 if (user->closure_tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL)); 107 PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets)); 108 PetscCall(DMGetCoordinateDim(dm, &cdim)); 109 PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim); 110 PetscCall(DMGetLocalSection(cdm, &s)); 111 PetscCall(VecGetArrayRead(X, &x)); 112 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element in %s order\n", cname, user->closure_tensor ? "tensor" : "bfs")); 113 for (PetscInt c = 0; c < Ncell; ++c) { 114 for (PetscInt v = 0; v < Ncl; ++v) { 115 PetscInt off = offsets[c * Ncl + v], dgdof; 116 const PetscScalar *vx = &x[off]; 117 118 if (isDG) { 119 PetscCall(PetscSectionGetDof(s, c, &dgdof)); 120 PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof); 121 } 122 switch (cdim) { 123 case 1: 124 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]))); 125 break; 126 case 2: 127 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]))); 128 break; 129 case 3: 130 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]))); 131 } 132 } 133 } 134 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout)); 135 PetscCall(VecRestoreArrayRead(X, &x)); 136 PetscCall(PetscFree(offsets)); 137 PetscCall(DMGetLocalToGlobalMapping(cdm, <og)); 138 PetscCall(ISLocalToGlobalMappingViewFromOptions(ltog, NULL, "-coord_ltog_view")); 139 } 140 PetscFunctionReturn(PETSC_SUCCESS); 141 } 142 143 int main(int argc, char **argv) 144 { 145 DM dm; 146 AppCtx user; 147 PetscInt depth; 148 149 PetscFunctionBeginUser; 150 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 151 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 152 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 153 PetscCall(SetupDiscretization(dm, &user)); 154 PetscCall(CheckOffsets(dm, &user, NULL, 0, 0)); 155 PetscCall(DMPlexGetDepth(dm, &depth)); 156 if (depth > 1) PetscCall(CheckOffsets(dm, &user, "Face Sets", user.check_face, 1)); 157 PetscCall(DMDestroy(&dm)); 158 PetscCall(PetscFinalize()); 159 return 0; 160 } 161 162 /*TEST 163 164 test: 165 suffix: 0 166 requires: triangle 167 args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view 168 169 test: 170 suffix: 1 171 args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \ 172 -dm_view -offsets_view 173 174 test: 175 suffix: cg_2d 176 args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \ 177 -dm_view -offsets_view 178 179 test: 180 suffix: 1d_sfc 181 args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_shape zbox -dm_plex_box_faces 3 1 -dm_view -coord_ltog_view 182 183 test: 184 suffix: 2d_sfc 185 nsize: 2 186 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 187 188 test: 189 suffix: 2d_sfc_periodic 190 nsize: 2 191 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 192 193 testset: 194 args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_shape zbox -dm_plex_box_faces 3,2 -petscspace_degree 1 -dm_plex_box_bd none,periodic -dm_view ::ascii_info_detail -closure_tensor 195 nsize: 2 196 test: 197 suffix: 2d_sfc_periodic_stranded 198 args: -dm_distribute 0 199 test: 200 suffix: 2d_sfc_periodic_stranded_dist 201 args: -dm_distribute 1 -petscpartitioner_type simple 202 203 TEST*/ 204