1 static char help[] = "Tests for creation of submeshes\n\n"; 2 3 #include <petscdmplex.h> 4 5 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm) 6 { 7 PetscFunctionBeginUser; 8 PetscCall(DMCreate(comm, dm)); 9 PetscCall(DMSetType(*dm, DMPLEX)); 10 PetscCall(DMSetFromOptions(*dm)); 11 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 12 PetscFunctionReturn(PETSC_SUCCESS); 13 } 14 15 // Label half of the cells 16 static PetscErrorCode CreateHalfCellsLabel(DM dm, PetscBool lower, DMLabel *label) 17 { 18 PetscInt cStart, cEnd, cStartSub, cEndSub; 19 20 PetscFunctionBeginUser; 21 PetscCall(DMCreateLabel(dm, "cells")); 22 PetscCall(DMGetLabel(dm, "cells", label)); 23 PetscCall(DMLabelClearStratum(*label, 1)); 24 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 25 if (lower) { 26 cStartSub = cStart; 27 cEndSub = cEnd / 2; 28 } else { 29 cStartSub = cEnd / 2; 30 cEndSub = cEnd; 31 } 32 for (PetscInt c = cStartSub; c < cEndSub; ++c) PetscCall(DMLabelSetValue(*label, c, 1)); 33 PetscCall(DMPlexLabelComplete(dm, *label)); 34 PetscFunctionReturn(PETSC_SUCCESS); 35 } 36 37 // Label everything on the right half of the domain 38 static PetscErrorCode CreateHalfDomainLabel(DM dm, PetscBool lower, PetscReal height, DMLabel *label) 39 { 40 PetscReal centroid[3]; 41 PetscInt cStart, cEnd, cdim; 42 43 PetscFunctionBeginUser; 44 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 45 PetscCall(DMCreateLabel(dm, "cells")); 46 PetscCall(DMGetLabel(dm, "cells", label)); 47 PetscCall(DMLabelClearStratum(*label, 1)); 48 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 49 PetscCall(DMGetCoordinateDim(dm, &cdim)); 50 for (PetscInt c = cStart; c < cEnd; ++c) { 51 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 52 if (height > 0.0 && PetscAbsReal(centroid[1] - height) > PETSC_SMALL) continue; 53 if (lower) { 54 if (centroid[0] < 0.5) PetscCall(DMLabelSetValue(*label, c, 1)); 55 } else { 56 if (centroid[0] > 0.5) PetscCall(DMLabelSetValue(*label, c, 1)); 57 } 58 } 59 PetscCall(DMPlexLabelComplete(dm, *label)); 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 63 // Create a line of faces at a given x value 64 static PetscErrorCode CreateLineLabel(DM dm, PetscReal x, DMLabel *label) 65 { 66 PetscReal centroid[3]; 67 PetscInt fStart, fEnd; 68 69 PetscFunctionBeginUser; 70 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 71 PetscCall(DMCreateLabel(dm, "faces")); 72 PetscCall(DMGetLabel(dm, "faces", label)); 73 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 74 for (PetscInt f = fStart; f < fEnd; ++f) { 75 PetscCall(DMPlexComputeCellGeometryFVM(dm, f, NULL, centroid, NULL)); 76 if (PetscAbsReal(centroid[0] - x) < PETSC_SMALL) PetscCall(DMLabelSetValue(*label, f, 1)); 77 } 78 PetscCall(DMPlexLabelComplete(dm, *label)); 79 PetscFunctionReturn(PETSC_SUCCESS); 80 } 81 82 static PetscErrorCode CreateVolumeSubmesh(DM dm, PetscBool domain, PetscBool lower, PetscReal height, DM *subdm) 83 { 84 DMLabel label, map; 85 86 PetscFunctionBegin; 87 if (domain) PetscCall(CreateHalfDomainLabel(dm, lower, height, &label)); 88 else PetscCall(CreateHalfCellsLabel(dm, lower, &label)); 89 PetscCall(DMPlexFilter(dm, label, 1, PETSC_FALSE, PETSC_FALSE, PetscObjectComm((PetscObject)dm), NULL, subdm)); 90 PetscCall(PetscObjectSetName((PetscObject)*subdm, "Submesh")); 91 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subdm, "sub_")); 92 PetscCall(DMViewFromOptions(*subdm, NULL, "-dm_view")); 93 PetscCall(DMPlexGetSubpointMap(*subdm, &map)); 94 PetscCall(PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view")); 95 PetscFunctionReturn(PETSC_SUCCESS); 96 } 97 98 static PetscErrorCode TestBoundaryField(DM dm) 99 { 100 DM subdm; 101 DMLabel label, map; 102 PetscFV fvm; 103 Vec gv; 104 PetscInt cdim; 105 106 PetscFunctionBeginUser; 107 PetscCall(CreateLineLabel(dm, 0.5, &label)); 108 PetscCall(DMPlexFilter(dm, label, 1, PETSC_FALSE, PETSC_FALSE, PetscObjectComm((PetscObject)dm), NULL, &subdm)); 109 PetscCall(PetscObjectSetName((PetscObject)subdm, "Submesh")); 110 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subdm, "sub_")); 111 PetscCall(DMViewFromOptions(subdm, NULL, "-dm_view")); 112 PetscCall(DMPlexGetSubpointMap(subdm, &map)); 113 PetscCall(PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view")); 114 115 PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)subdm), &fvm)); 116 PetscCall(PetscObjectSetName((PetscObject)fvm, "testField")); 117 PetscCall(PetscFVSetNumComponents(fvm, 1)); 118 PetscCall(DMGetCoordinateDim(subdm, &cdim)); 119 PetscCall(PetscFVSetSpatialDimension(fvm, cdim)); 120 PetscCall(PetscFVSetFromOptions(fvm)); 121 122 PetscCall(DMAddField(subdm, NULL, (PetscObject)fvm)); 123 PetscCall(PetscFVDestroy(&fvm)); 124 PetscCall(DMCreateDS(subdm)); 125 126 PetscCall(DMCreateGlobalVector(subdm, &gv)); 127 PetscCall(PetscObjectSetName((PetscObject)gv, "potential")); 128 PetscCall(VecSet(gv, 1.)); 129 PetscCall(VecViewFromOptions(gv, NULL, "-vec_view")); 130 PetscCall(VecDestroy(&gv)); 131 PetscCall(DMDestroy(&subdm)); 132 PetscFunctionReturn(PETSC_SUCCESS); 133 } 134 135 int main(int argc, char **argv) 136 { 137 DM dm, subdm; 138 PetscReal height = -1.0; 139 PetscBool volume = PETSC_TRUE, domain = PETSC_FALSE; 140 141 PetscFunctionBeginUser; 142 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 143 PetscCall(PetscOptionsGetReal(NULL, NULL, "-height", &height, NULL)); 144 PetscCall(PetscOptionsGetBool(NULL, NULL, "-volume", &volume, NULL)); 145 PetscCall(PetscOptionsGetBool(NULL, NULL, "-domain", &domain, NULL)); 146 147 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm)); 148 if (volume) { 149 PetscCall(CreateVolumeSubmesh(dm, domain, PETSC_TRUE, height, &subdm)); 150 PetscCall(DMSetFromOptions(subdm)); 151 PetscCall(DMDestroy(&subdm)); 152 PetscCall(CreateVolumeSubmesh(dm, domain, PETSC_FALSE, height, &subdm)); 153 PetscCall(DMSetFromOptions(subdm)); 154 PetscCall(DMDestroy(&subdm)); 155 } else { 156 PetscCall(TestBoundaryField(dm)); 157 } 158 PetscCall(DMDestroy(&dm)); 159 PetscCall(PetscFinalize()); 160 return 0; 161 } 162 163 /*TEST 164 165 test: 166 suffix: 0 167 requires: triangle 168 args: -dm_coord_space 0 -sub_dm_plex_check_all \ 169 -dm_view ascii::ascii_info_detail -sub_dm_view ascii::ascii_info_detail -map_view 170 171 test: 172 suffix: 0_vtk 173 requires: triangle 174 args: -dm_coord_space 0 -sub_dm_plex_check_all \ 175 -dm_view vtk: -sub_dm_view vtk: -map_view 176 177 # These tests check that filtering is stable when boundary point ownership could change, so it needs 3 processes 178 testset: 179 nsize: 3 180 requires: parmetis 181 args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -petscpartitioner_type parmetis -dm_distribute_overlap 1 -sub_dm_distribute 0 \ 182 -sub_dm_plex_check_all -dm_view -sub_dm_view 183 184 test: 185 suffix: 1 186 187 test: 188 suffix: 2 189 args: -domain 190 191 # This set tests that global numberings can be made when some strata are missing on a process 192 testset: 193 nsize: 3 194 requires: hdf5 195 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \ 196 -sub_dm_plex_check_all -sub_dm_view hdf5:subdm.h5 197 198 test: 199 suffix: 3 200 args: -domain -height 0.625 201 output_file: output/empty.out 202 203 # This set tests that global numberings can be made when some strata are missing on a process 204 testset: 205 nsize: 3 206 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,4 -petscpartitioner_type simple -sub_dm_distribute 0 \ 207 -sub_dm_plex_check_all -sub_dm_view {{vtk:subdm.vtu :subdm.txt :subdm_d.txt:ascii_info_detail}} 208 209 test: 210 suffix: 3_vtk 211 args: -domain -height 0.625 212 output_file: output/empty.out 213 214 # This test checks whether filter can extract a lower-dimensional manifold and output a field on it 215 testset: 216 requires: hdf5 217 args: -volume 0 -dm_plex_simplex 0 -sub_dm_view hdf5:subdm.h5 -vec_view hdf5:subdm.h5::append 218 output_file: output/empty.out 219 220 test: 221 suffix: surface_2d 222 args: -dm_plex_box_faces 5,5 223 224 test: 225 suffix: surface_3d 226 args: -dm_plex_box_faces 5,5,5 227 228 # This test checks that if cells are present in the SF, the dm is marked as having an overlap 229 test: 230 suffix: surface_2d_2 231 nsize: 3 232 args: -dm_plex_box_faces 5,5 -domain -height 0.625 233 234 TEST*/ 235