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(0); 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) {cStartSub = cStart; cEndSub = cEnd/2;} 26 else {cStartSub = cEnd/2; cEndSub = cEnd;} 27 for (PetscInt c = cStartSub; c < cEndSub; ++c) PetscCall(DMLabelSetValue(*label, c, 1)); 28 PetscCall(DMPlexLabelComplete(dm, *label)); 29 PetscFunctionReturn(0); 30 } 31 32 // Label everything on the right half of the domain 33 static PetscErrorCode CreateHalfDomainLabel(DM dm, PetscBool lower, DMLabel *label) 34 { 35 PetscReal centroid[3]; 36 PetscInt cStart, cEnd, cdim; 37 38 PetscFunctionBeginUser; 39 PetscCall(DMCreateLabel(dm, "cells")); 40 PetscCall(DMGetLabel(dm, "cells", label)); 41 PetscCall(DMLabelClearStratum(*label, 1)); 42 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 43 PetscCall(DMGetCoordinateDim(dm, &cdim)); 44 for (PetscInt c = cStart; c < cEnd; ++c) { 45 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL)); 46 if (lower) {if (centroid[0] < 0.5) PetscCall(DMLabelSetValue(*label, c, 1));} 47 else {if (centroid[0] > 0.5) PetscCall(DMLabelSetValue(*label, c, 1));} 48 } 49 PetscCall(DMPlexLabelComplete(dm, *label)); 50 PetscFunctionReturn(0); 51 } 52 53 PetscErrorCode CreateSubmesh(DM dm, PetscBool domain, PetscBool lower, DM *subdm) 54 { 55 DMLabel label, map; 56 57 PetscFunctionBegin; 58 if (domain) PetscCall(CreateHalfDomainLabel(dm, lower, &label)); 59 else PetscCall(CreateHalfCellsLabel(dm, lower, &label)); 60 PetscCall(DMPlexFilter(dm, label, 1, subdm)); 61 PetscCall(PetscObjectSetName((PetscObject) *subdm, "Submesh")); 62 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *subdm, "sub_")); 63 PetscCall(DMViewFromOptions(*subdm, NULL, "-dm_view")); 64 PetscCall(DMPlexGetSubpointMap(*subdm, &map)); 65 PetscCall(PetscObjectViewFromOptions((PetscObject) map, NULL, "-map_view")); 66 PetscFunctionReturn(0); 67 } 68 69 int main(int argc, char **argv) 70 { 71 DM dm, subdm; 72 PetscBool domain = PETSC_FALSE; 73 74 PetscFunctionBeginUser; 75 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 76 PetscCall(PetscOptionsGetBool(NULL, NULL, "-domain", &domain, NULL)); 77 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm)); 78 PetscCall(CreateSubmesh(dm, domain, PETSC_TRUE, &subdm)); 79 PetscCall(DMSetFromOptions(subdm)); 80 PetscCall(DMDestroy(&subdm)); 81 PetscCall(CreateSubmesh(dm, domain, PETSC_FALSE, &subdm)); 82 PetscCall(DMSetFromOptions(subdm)); 83 PetscCall(DMDestroy(&subdm)); 84 PetscCall(DMDestroy(&dm)); 85 PetscCall(PetscFinalize()); 86 return 0; 87 } 88 89 /*TEST 90 91 test: 92 suffix: 0 93 requires: triangle 94 args: -dm_coord_space 0 -sub_dm_plex_check_all \ 95 -dm_view ascii::ascii_info_detail -sub_dm_view ascii::ascii_info_detail -map_view 96 97 # These tests check that filtering is stable when boundary point ownership could change, so it needs 3 processes 98 testset: 99 nsize: 3 100 requires: parmetis 101 args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -petscpartitioner_type parmetis -dm_distribute_overlap 1 -sub_dm_distribute 0 \ 102 -sub_dm_plex_check_all -dm_view -sub_dm_view 103 104 test: 105 suffix: 1 106 107 test: 108 suffix: 2 109 args: -domain 110 111 TEST*/ 112