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 PetscErrorCode CreateSubmesh(DM dm, PetscBool domain, PetscBool lower, DM *subdm) { 59 DMLabel label, map; 60 61 PetscFunctionBegin; 62 if (domain) PetscCall(CreateHalfDomainLabel(dm, lower, &label)); 63 else PetscCall(CreateHalfCellsLabel(dm, lower, &label)); 64 PetscCall(DMPlexFilter(dm, label, 1, subdm)); 65 PetscCall(PetscObjectSetName((PetscObject)*subdm, "Submesh")); 66 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subdm, "sub_")); 67 PetscCall(DMViewFromOptions(*subdm, NULL, "-dm_view")); 68 PetscCall(DMPlexGetSubpointMap(*subdm, &map)); 69 PetscCall(PetscObjectViewFromOptions((PetscObject)map, NULL, "-map_view")); 70 PetscFunctionReturn(0); 71 } 72 73 int main(int argc, char **argv) { 74 DM dm, subdm; 75 PetscBool domain = PETSC_FALSE; 76 77 PetscFunctionBeginUser; 78 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 79 PetscCall(PetscOptionsGetBool(NULL, NULL, "-domain", &domain, NULL)); 80 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm)); 81 PetscCall(CreateSubmesh(dm, domain, PETSC_TRUE, &subdm)); 82 PetscCall(DMSetFromOptions(subdm)); 83 PetscCall(DMDestroy(&subdm)); 84 PetscCall(CreateSubmesh(dm, domain, PETSC_FALSE, &subdm)); 85 PetscCall(DMSetFromOptions(subdm)); 86 PetscCall(DMDestroy(&subdm)); 87 PetscCall(DMDestroy(&dm)); 88 PetscCall(PetscFinalize()); 89 return 0; 90 } 91 92 /*TEST 93 94 test: 95 suffix: 0 96 requires: triangle 97 args: -dm_coord_space 0 -sub_dm_plex_check_all \ 98 -dm_view ascii::ascii_info_detail -sub_dm_view ascii::ascii_info_detail -map_view 99 100 # These tests check that filtering is stable when boundary point ownership could change, so it needs 3 processes 101 testset: 102 nsize: 3 103 requires: parmetis 104 args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -petscpartitioner_type parmetis -dm_distribute_overlap 1 -sub_dm_distribute 0 \ 105 -sub_dm_plex_check_all -dm_view -sub_dm_view 106 107 test: 108 suffix: 1 109 110 test: 111 suffix: 2 112 args: -domain 113 114 TEST*/ 115