1 static char help[] = "Tests interpolation and output of hybrid meshes\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petscsf.h> 5 6 /* Much can be learned using 7 -rd_dm_view -rd2_dm_view -rd_section_view -rd_vec_view -rd2_section_view */ 8 9 static PetscErrorCode redistribute_vec(DM dist_dm, PetscSF sf, Vec *v) { 10 DM dm, dist_v_dm; 11 PetscSection section, dist_section; 12 Vec dist_v; 13 PetscMPIInt rank, size, p; 14 15 PetscFunctionBegin; 16 PetscCall(VecGetDM(*v, &dm)); 17 PetscCall(DMGetLocalSection(dm, §ion)); 18 PetscCall(DMViewFromOptions(dm, NULL, "-rd_dm_view")); 19 PetscCall(DMViewFromOptions(dist_dm, NULL, "-rd2_dm_view")); 20 21 PetscCall(DMClone(dm, &dist_v_dm)); 22 PetscCall(VecCreate(PetscObjectComm((PetscObject)*v), &dist_v)); 23 PetscCall(VecSetDM(dist_v, dist_v_dm)); 24 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)*v), &dist_section)); 25 PetscCall(DMSetLocalSection(dist_v_dm, dist_section)); 26 27 PetscCall(PetscObjectViewFromOptions((PetscObject)section, NULL, "-rd_section_view")); 28 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 29 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size)); 30 for (p = 0; p < size; ++p) { 31 if (p == rank) PetscCall(PetscObjectViewFromOptions((PetscObject)*v, NULL, "-rd_vec_view")); 32 PetscCall(PetscBarrier((PetscObject)dm)); 33 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 34 } 35 PetscCall(DMPlexDistributeField(dm, sf, section, *v, dist_section, dist_v)); 36 for (p = 0; p < size; ++p) { 37 if (p == rank) { 38 PetscCall(PetscObjectViewFromOptions((PetscObject)dist_section, NULL, "-rd2_section_view")); 39 PetscCall(PetscObjectViewFromOptions((PetscObject)dist_v, NULL, "-rd2_vec_view")); 40 } 41 PetscCall(PetscBarrier((PetscObject)dm)); 42 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 43 } 44 45 PetscCall(PetscSectionDestroy(&dist_section)); 46 PetscCall(DMDestroy(&dist_v_dm)); 47 48 PetscCall(VecDestroy(v)); 49 *v = dist_v; 50 PetscFunctionReturn(0); 51 } 52 53 static PetscErrorCode dm_view_geometry(DM dm, Vec cell_geom, Vec face_geom) { 54 DM cell_dm, face_dm; 55 PetscSection cell_section, face_section; 56 const PetscScalar *cell_array, *face_array; 57 const PetscInt *cells; 58 PetscInt c, start_cell, end_cell; 59 PetscInt f, start_face, end_face; 60 PetscInt supportSize, offset; 61 PetscMPIInt rank; 62 63 PetscFunctionBegin; 64 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 65 66 /* cells */ 67 PetscCall(DMPlexGetHeightStratum(dm, 0, &start_cell, &end_cell)); 68 PetscCall(VecGetDM(cell_geom, &cell_dm)); 69 PetscCall(DMGetLocalSection(cell_dm, &cell_section)); 70 PetscCall(VecGetArrayRead(cell_geom, &cell_array)); 71 72 for (c = start_cell; c < end_cell; ++c) { 73 const PetscFVCellGeom *geom; 74 PetscCall(PetscSectionGetOffset(cell_section, c, &offset)); 75 geom = (PetscFVCellGeom *)&cell_array[offset]; 76 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "rank %d c %" PetscInt_FMT " centroid %g,%g,%g vol %g\n", rank, c, (double)geom->centroid[0], (double)geom->centroid[1], (double)geom->centroid[2], (double)geom->volume)); 77 } 78 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL)); 79 PetscCall(VecRestoreArrayRead(cell_geom, &cell_array)); 80 81 /* faces */ 82 PetscCall(DMPlexGetHeightStratum(dm, 1, &start_face, &end_face)); 83 PetscCall(VecGetDM(face_geom, &face_dm)); 84 PetscCall(DMGetLocalSection(face_dm, &face_section)); 85 PetscCall(VecGetArrayRead(face_geom, &face_array)); 86 for (f = start_face; f < end_face; ++f) { 87 PetscCall(DMPlexGetSupport(dm, f, &cells)); 88 PetscCall(DMPlexGetSupportSize(dm, f, &supportSize)); 89 if (supportSize > 1) { 90 PetscCall(PetscSectionGetOffset(face_section, f, &offset)); 91 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "rank %d f %" PetscInt_FMT " cells %" PetscInt_FMT ",%" PetscInt_FMT " normal %g,%g,%g centroid %g,%g,%g\n", rank, f, cells[0], cells[1], (double)PetscRealPart(face_array[offset + 0]), (double)PetscRealPart(face_array[offset + 1]), (double)PetscRealPart(face_array[offset + 2]), (double)PetscRealPart(face_array[offset + 3]), (double)PetscRealPart(face_array[offset + 4]), (double)PetscRealPart(face_array[offset + 5]))); 92 } 93 } 94 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL)); 95 PetscCall(VecRestoreArrayRead(cell_geom, &cell_array)); 96 PetscFunctionReturn(0); 97 } 98 99 int main(int argc, char **argv) { 100 DM dm, redist_dm; 101 PetscPartitioner part; 102 PetscSF redist_sf; 103 Vec cell_geom, face_geom; 104 PetscInt overlap2 = 1; 105 106 PetscFunctionBeginUser; 107 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 108 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 109 PetscCall(DMSetType(dm, DMPLEX)); 110 PetscCall(DMSetFromOptions(dm)); 111 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 112 113 PetscCall(DMPlexComputeGeometryFVM(dm, &cell_geom, &face_geom)); 114 PetscCall(dm_view_geometry(dm, cell_geom, face_geom)); 115 116 /* redistribute */ 117 PetscCall(DMPlexGetPartitioner(dm, &part)); 118 PetscCall(PetscPartitionerSetFromOptions(part)); 119 PetscCall(PetscOptionsGetInt(NULL, NULL, "-overlap2", &overlap2, NULL)); 120 PetscCall(DMPlexDistribute(dm, overlap2, &redist_sf, &redist_dm)); 121 if (redist_dm) { 122 PetscCall(redistribute_vec(redist_dm, redist_sf, &cell_geom)); 123 PetscCall(redistribute_vec(redist_dm, redist_sf, &face_geom)); 124 PetscCall(PetscObjectViewFromOptions((PetscObject)redist_sf, NULL, "-rd2_sf_view")); 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "redistributed:\n")); 126 PetscCall(dm_view_geometry(redist_dm, cell_geom, face_geom)); 127 } 128 129 PetscCall(VecDestroy(&cell_geom)); 130 PetscCall(VecDestroy(&face_geom)); 131 PetscCall(PetscSFDestroy(&redist_sf)); 132 PetscCall(DMDestroy(&redist_dm)); 133 PetscCall(DMDestroy(&dm)); 134 PetscCall(PetscFinalize()); 135 return 0; 136 } 137 138 /*TEST 139 140 test: 141 suffix: 0 142 nsize: 3 143 args: -dm_plex_dim 3 -dm_plex_box_faces 8,1,1 -dm_plex_simplex 0 -dm_plex_adj_cone 1 -dm_plex_adj_closure 0 -petscpartitioner_type simple -dm_distribute_overlap 1 -overlap2 1 144 145 TEST*/ 146