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