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