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