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