1 static char help[] = "Element closure restrictions in tensor/lexicographic/spectral-element ordering using DMPlex\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petscds.h> 5 6 static PetscErrorCode ViewOffsets(DM dm, Vec X) 7 { 8 PetscInt num_elem, elem_size, num_comp, num_dof; 9 PetscInt *elem_restr_offsets; 10 const PetscScalar *x = NULL; 11 const char *name; 12 13 PetscFunctionBegin; 14 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 15 PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets)); 16 PetscCall(PetscPrintf(PETSC_COMM_SELF, "DM %s offsets: num_elem %" PetscInt_FMT ", size %" PetscInt_FMT ", comp %" PetscInt_FMT ", dof %" PetscInt_FMT "\n", name, num_elem, elem_size, num_comp, num_dof)); 17 if (X) PetscCall(VecGetArrayRead(X, &x)); 18 for (PetscInt c = 0; c < num_elem; c++) { 19 PetscCall(PetscIntView(elem_size, &elem_restr_offsets[c * elem_size], PETSC_VIEWER_STDOUT_SELF)); 20 if (x) { 21 for (PetscInt i = 0; i < elem_size; i++) PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c * elem_size + i]], PETSC_VIEWER_STDERR_SELF)); 22 } 23 } 24 if (X) PetscCall(VecRestoreArrayRead(X, &x)); 25 PetscCall(PetscFree(elem_restr_offsets)); 26 PetscFunctionReturn(PETSC_SUCCESS); 27 } 28 29 int main(int argc, char **argv) 30 { 31 DM dm; 32 PetscSection section; 33 PetscFE fe; 34 PetscInt dim, Nf = 1, c, cStart, cEnd; 35 PetscBool view_coord = PETSC_FALSE, tensor = PETSC_TRUE, project = PETSC_FALSE; 36 37 PetscFunctionBeginUser; 38 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 39 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Tensor closure restrictions", "DMPLEX"); 40 PetscCall(PetscOptionsBool("-closure_tensor", "Apply DMPlexSetClosurePermutationTensor", "ex8.c", tensor, &tensor, NULL)); 41 PetscCall(PetscOptionsBool("-project_coordinates", "Call DMSetCoordinateDisc() explicitly", "ex8.c", project, &project, NULL)); 42 PetscCall(PetscOptionsBool("-view_coord", "View coordinates of element closures", "ex8.c", view_coord, &view_coord, NULL)); 43 PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields to use", "ex8.c", Nf, &Nf, NULL, 1)); 44 PetscOptionsEnd(); 45 46 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 47 PetscCall(DMSetType(dm, DMPLEX)); 48 PetscCall(DMSetFromOptions(dm)); 49 if (project) { 50 PetscFE fe_coords; 51 PetscInt cdim; 52 PetscCall(DMGetCoordinateDim(dm, &cdim)); 53 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, cdim, cdim, PETSC_FALSE, 1, 1, &fe_coords)); 54 PetscCall(DMSetCoordinateDisc(dm, fe_coords, PETSC_FALSE, PETSC_TRUE)); 55 PetscCall(PetscFEDestroy(&fe_coords)); 56 } 57 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 58 PetscCall(DMGetDimension(dm, &dim)); 59 60 if (Nf == 1) { 61 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe)); 62 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 63 PetscCall(PetscFEDestroy(&fe)); 64 } else { 65 for (PetscInt f = 0; f < Nf; ++f) { 66 char prefix[16]; 67 PetscCall(PetscSNPrintf(prefix, 16, "f%" PetscInt_FMT "_", f)); 68 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, prefix, PETSC_DETERMINE, &fe)); 69 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 70 PetscCall(PetscFEDestroy(&fe)); 71 } 72 } 73 PetscCall(DMCreateDS(dm)); 74 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 75 PetscCall(DMGetLocalSection(dm, §ion)); 76 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 77 for (c = cStart; c < cEnd; c++) { 78 PetscInt numindices, *indices; 79 PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL)); 80 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "\n", c - cStart)); 81 PetscCall(PetscIntView(numindices, indices, PETSC_VIEWER_STDOUT_SELF)); 82 PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &numindices, &indices, NULL, NULL)); 83 } 84 { 85 DMLabel domain_label; 86 IS valueIS; 87 const PetscInt *values; 88 PetscInt Nv; 89 90 PetscCall(DMGetLabel(dm, "Face Sets", &domain_label)); 91 PetscCall(DMLabelGetNumValues(domain_label, &Nv)); 92 // Check that discretization bas any values on faces 93 { 94 PetscDS ds; 95 PetscFE fe; 96 IS fieldIS; 97 const PetscInt *fields; 98 PetscInt Nf, dmf = 0, dsf = -1; 99 100 PetscCall(DMGetRegionDS(dm, domain_label, &fieldIS, &ds, NULL)); 101 // Translate DM field number to DS field number 102 PetscCall(ISGetIndices(fieldIS, &fields)); 103 PetscCall(ISGetSize(fieldIS, &Nf)); 104 for (PetscInt f = 0; f < Nf; ++f) { 105 if (dmf == fields[f]) { 106 dsf = f; 107 break; 108 } 109 } 110 PetscCall(ISRestoreIndices(fieldIS, &fields)); 111 PetscCall(PetscDSGetDiscretization(ds, dsf, (PetscObject *)&fe)); 112 PetscCall(PetscFEGetHeightSubspace(fe, 1, &fe)); 113 if (!fe) Nv = 0; 114 } 115 PetscCall(DMLabelGetValueIS(domain_label, &valueIS)); 116 PetscCall(ISGetIndices(valueIS, &values)); 117 for (PetscInt v = 0; v < Nv; ++v) { 118 PetscInt *elem_restr_offsets_face; 119 PetscInt num_elem, elem_size, num_comp, num_dof; 120 121 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, values[v], 1, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets_face)); 122 PetscCall(PetscPrintf(PETSC_COMM_SELF, "========= Face restriction; marker: %" PetscInt_FMT " ========\n", values[v])); 123 PetscCall(PetscPrintf(PETSC_COMM_SELF, "num_elem: %" PetscInt_FMT ", elem_size: %" PetscInt_FMT ", num_dof: %" PetscInt_FMT "\n", num_elem, elem_size, num_dof)); 124 for (PetscInt c = 0; c < num_elem; c++) PetscCall(PetscIntView(elem_size, &elem_restr_offsets_face[c * elem_size], PETSC_VIEWER_STDOUT_SELF)); 125 PetscCall(PetscFree(elem_restr_offsets_face)); 126 } 127 PetscCall(ISRestoreIndices(valueIS, &values)); 128 PetscCall(ISDestroy(&valueIS)); 129 } 130 if (view_coord) { 131 DM cdm, cell_dm; 132 Vec X; 133 PetscInt cdim; 134 PetscBool sparseLocalize; 135 136 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 137 PetscCall(DMGetCoordinateDim(dm, &cdim)); 138 PetscCall(DMGetCoordinateDM(dm, &cdm)); 139 PetscCall(PetscObjectSetName((PetscObject)cdm, "coords")); 140 if (tensor) PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL)); 141 PetscCall(DMLocalizeCoordinates(dm)); 142 PetscCall(DMGetSparseLocalize(dm, &sparseLocalize)); 143 for (c = cStart; c < cEnd; ++c) { 144 const PetscScalar *array; 145 PetscScalar *x = NULL; 146 PetscInt ndof; 147 PetscBool isDG; 148 149 PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 150 PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof not divisible by cdim"); 151 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT " coordinates\n", c - cStart)); 152 for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim, &x[i], PETSC_VIEWER_STDOUT_SELF)); 153 PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x)); 154 } 155 PetscCall(ViewOffsets(dm, NULL)); 156 PetscCall(DMGetCoordinatesLocal(dm, &X)); 157 PetscCall(ViewOffsets(cdm, X)); 158 PetscCall(DMGetCellCoordinateDM(dm, &cell_dm)); 159 PetscCall(PetscObjectSetName((PetscObject)cell_dm, "cell coords")); 160 if (cell_dm && !sparseLocalize) { 161 PetscCall(DMGetCellCoordinatesLocal(dm, &X)); 162 PetscCall(ViewOffsets(cell_dm, X)); 163 } 164 } 165 PetscCall(DMDestroy(&dm)); 166 PetscCall(PetscFinalize()); 167 return 0; 168 } 169 170 /*TEST 171 172 test: 173 suffix: 1d_q2 174 args: -dm_plex_dim 1 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2 175 test: 176 suffix: 2d_q1 177 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 178 test: 179 suffix: 2d_q1d 180 args: -dm_plex_dim 2 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 -petscdualspace_lagrange_continuity 0 181 test: 182 suffix: 2d_q1_q1d 183 args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 1 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0 184 test: 185 suffix: 2d_q2 186 args: -dm_plex_dim 2 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2 187 test: 188 suffix: 2d_q2_q1 189 args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 190 test: 191 suffix: 2d_q2_p0 192 args: -dm_plex_dim 2 -num_fields 2 -f0_petscspace_degree 2 -f1_petscspace_degree 0 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -f1_petscdualspace_lagrange_continuity 0 193 test: 194 suffix: 2d_q2_q1d 195 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \ 196 -f0_petscspace_degree 2 \ 197 -f1_petscspace_degree 1 -f1_petscdualspace_lagrange_continuity 0 198 test: 199 suffix: 2d_q2_p1d 200 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,1 -num_fields 2 \ 201 -f0_petscspace_degree 2 \ 202 -f1_petscspace_degree 1 -f1_petscspace_poly_tensor 0 -f1_petscdualspace_lagrange_continuity 0 203 test: 204 suffix: 2d_q3 205 args: -dm_plex_dim 2 -petscspace_degree 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 206 207 test: 208 suffix: 3d_q1 209 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 210 test: 211 suffix: 1d_q1_periodic 212 requires: !complex 213 args: -dm_plex_dim 1 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3 -dm_plex_box_bd periodic -dm_view -view_coord 214 testset: 215 suffix: 2d_q1_periodic 216 requires: !complex 217 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 3,2 -dm_plex_box_bd periodic,none \ 218 -petscspace_degree 1 -dm_view -view_coord 219 220 test: 221 test: 222 suffix: sparse 223 args: -dm_sparse_localize false -dm_localize 0 224 225 test: 226 suffix: 3d_q1_periodic 227 requires: !complex 228 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 3,2,1 -dm_plex_box_bd periodic,none,none -dm_view -view_coord 229 test: 230 suffix: 3d_q1_periodic_project 231 requires: !complex 232 args: -dm_plex_dim 3 -petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,3 -dm_plex_box_bd none,none,periodic -dm_view -view_coord -project_coordinates 233 234 test: 235 suffix: 3d_q2_periodic # not actually periodic because only 2 cells 236 args: -dm_plex_dim 3 -petscspace_degree 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2 -dm_plex_box_bd periodic,none,periodic -dm_view 237 238 TEST*/ 239