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)); 88 cname = "Coordinates"; 89 PetscCall(DMGetCoordinatesLocal(dm, &X)); 90 } else { 91 isDG = PETSC_TRUE; 92 cname = "DG Coordinates"; 93 PetscCall(DMGetCellCoordinatesLocal(dm, &X)); 94 } 95 if (isDG && height) PetscFunctionReturn(0); 96 if (domain_name) PetscCall(DMGetLabel(cdm, domain_name, &domain_label)); 97 PetscCall(DMPlexGetLocalOffsets(cdm, domain_label, label_value, height, 0, &Ncell, &Ncl, &Nc, &n, &offsets)); 98 PetscCall(DMGetCoordinateDim(dm, &cdim)); 99 PetscCheck(Nc == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometric dimension %" PetscInt_FMT " should be %" PetscInt_FMT, Nc, cdim); 100 PetscCall(DMGetLocalSection(cdm, &s)); 101 PetscCall(VecGetArrayRead(X, &x)); 102 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s by element\n", cname)); 103 for (PetscInt c = 0; c < Ncell; ++c) { 104 for (PetscInt v = 0; v < Ncl; ++v) { 105 PetscInt off = offsets[c * Ncl + v], dgdof; 106 const PetscScalar *vx = &x[off]; 107 108 if (isDG) { 109 PetscCall(PetscSectionGetDof(s, c, &dgdof)); 110 PetscCheck(Ncl * Nc == dgdof, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Offset size %" PetscInt_FMT " should be %" PetscInt_FMT, Ncl * Nc, dgdof); 111 } 112 switch (cdim) { 113 case 2: 114 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]))); 115 break; 116 case 3: 117 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]))); 118 break; 119 } 120 } 121 } 122 PetscCall(VecRestoreArrayRead(X, &x)); 123 PetscCall(PetscFree(offsets)); 124 } 125 PetscFunctionReturn(0); 126 } 127 128 int main(int argc, char **argv) 129 { 130 DM dm; 131 AppCtx user; 132 133 PetscFunctionBeginUser; 134 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 135 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 136 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 137 PetscCall(SetupDiscretization(dm, &user)); 138 PetscCall(CheckOffsets(dm, NULL, 0, 0)); 139 PetscCall(CheckOffsets(dm, "Face Sets", 1, 1)); 140 PetscCall(DMDestroy(&dm)); 141 PetscCall(PetscFinalize()); 142 return 0; 143 } 144 145 /*TEST 146 147 test: 148 suffix: 0 149 requires: triangle 150 args: -dm_refine 1 -petscspace_degree 1 -dm_view -offsets_view 151 152 test: 153 suffix: 1 154 args: -dm_plex_simplex 0 -dm_plex_box_bd periodic,none -dm_plex_box_faces 3,3 -dm_sparse_localize 0 -petscspace_degree 1 \ 155 -dm_view -offsets_view 156 157 test: 158 suffix: cg_2d 159 args: -dm_plex_simplex 0 -dm_plex_box_bd none,none -dm_plex_box_faces 3,3 -petscspace_degree 1 \ 160 -dm_view -offsets_view 161 TEST*/ 162