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