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