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 61 PetscFunctionBeginUser; 62 if (domain_name) PetscCall(DMGetLabel(dm, domain_name, &domain_label)); 63 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)" : "")); 64 if (domain_name && !domain_label) PetscFunctionReturn(0); 65 // Offsets for cell closures 66 PetscCall(DMGetNumFields(dm, &Nf)); 67 for (f = 0; f < Nf; ++f) { 68 char name[PETSC_MAX_PATH_LEN]; 69 70 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, label_value, height, f, &Ncell, &Ncl, &Nc, &n, &offsets)); 71 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncell * Ncl, offsets, PETSC_OWN_POINTER, &offIS)); 72 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Field %" PetscInt_FMT " Offsets", f)); 73 PetscCall(PetscObjectSetName((PetscObject)offIS, name)); 74 PetscCall(ISViewFromOptions(offIS, NULL, "-offsets_view")); 75 PetscCall(ISDestroy(&offIS)); 76 } 77 // Offsets for coordinates 78 { 79 Vec X; 80 PetscSection s; 81 const PetscScalar *x; 82 const char *cname; 83 PetscInt cdim; 84 PetscBool isDG = PETSC_FALSE; 85 86 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 87 if (!cdm) { 88 PetscCall(DMGetCoordinateDM(dm, &cdm)); 89 cname = "Coordinates"; 90 PetscCall(DMGetCoordinatesLocal(dm, &X)); 91 } else { 92 isDG = PETSC_TRUE; 93 cname = "DG Coordinates"; 94 PetscCall(DMGetCellCoordinatesLocal(dm, &X)); 95 } 96 if (isDG && height) PetscFunctionReturn(0); 97 if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label)); 98 PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets)); 99 PetscCall(DMGetCoordinateDim(dm, &cdim)); 100 PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim); 101 PetscCall(DMGetLocalSection(cdm, &s)); 102 PetscCall(VecGetArrayRead(X, &x)); 103 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element\n", cname)); 104 for (PetscInt c = 0; c < Ncell; ++c) { 105 for (PetscInt v = 0; v < Ncl; ++v) { 106 PetscInt off = offsets[c * Ncl + v], dgdof; 107 const PetscScalar *vx = &x[off]; 108 109 if (isDG) { 110 PetscCall(PetscSectionGetDof(s, c, &dgdof)); 111 PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof); 112 } 113 switch (cdim) { 114 case 1: 115 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT "] %" PetscInt_FMT " <-- %2" PetscInt_FMT " (% 4.2f)\n", c, v, off, (double)PetscRealPart(vx[0]))); 116 break; 117 case 2: 118 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]))); 119 break; 120 case 3: 121 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]))); 122 } 123 } 124 } 125 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout)); 126 PetscCall(VecRestoreArrayRead(X, &x)); 127 PetscCall(PetscFree(offsets)); 128 } 129 PetscFunctionReturn(0); 130 } 131 132 int main(int argc, char **argv) 133 { 134 DM dm; 135 AppCtx user; 136 PetscInt depth; 137 138 PetscFunctionBeginUser; 139 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 140 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 141 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 142 PetscCall(SetupDiscretization(dm, &user)); 143 PetscCall(CheckOffsets(dm, NULL, 0, 0)); 144 PetscCall(DMPlexGetDepth(dm, &depth)); 145 if (depth > 1) PetscCall(CheckOffsets(dm, "Face Sets", 1, 1)); 146 PetscCall(DMDestroy(&dm)); 147 PetscCall(PetscFinalize()); 148 return 0; 149 } 150 151 /*TEST 152 153 test: 154 suffix: 0 155 requires: triangle 156 args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view 157 158 test: 159 suffix: 1 160 args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \ 161 -dm_view -offsets_view 162 163 test: 164 suffix: cg_2d 165 args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \ 166 -dm_view -offsets_view 167 168 test: 169 suffix: 1d_sfc 170 args: -dm_plex_simplex 0 -dm_plex_dim 1 -dm_plex_box_faces 3 -dm_plex_box_sfc 1 -dm_view 171 172 test: 173 suffix: 2d_sfc 174 nsize: 2 175 args: -dm_plex_simplex 0 -dm_plex_dim 2 -dm_plex_box_faces 4,3 -dm_plex_box_sfc -dm_distribute 0 -dm_view 176 177 TEST*/ 178