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 CHKERRQ(VecGetDM(*v, &dm)); 18 CHKERRQ(DMGetLocalSection(dm, §ion)); 19 CHKERRQ(DMViewFromOptions(dm, NULL, "-rd_dm_view")); 20 CHKERRQ(DMViewFromOptions(dist_dm, NULL, "-rd2_dm_view")); 21 22 CHKERRQ(DMClone(dm, &dist_v_dm)); 23 CHKERRQ(VecCreate(PetscObjectComm((PetscObject) *v), &dist_v)); 24 CHKERRQ(VecSetDM(dist_v, dist_v_dm)); 25 CHKERRQ(PetscSectionCreate(PetscObjectComm((PetscObject) *v), &dist_section)); 26 CHKERRQ(DMSetLocalSection(dist_v_dm, dist_section)); 27 28 CHKERRQ(PetscObjectViewFromOptions((PetscObject) section, NULL, "-rd_section_view")); 29 CHKERRMPI(MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank)); 30 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 31 for (p = 0; p < size; ++p) { 32 if (p == rank) { 33 CHKERRQ(PetscObjectViewFromOptions((PetscObject) *v, NULL, "-rd_vec_view"));} 34 CHKERRQ(PetscBarrier((PetscObject) dm)); 35 CHKERRQ(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 36 } 37 CHKERRQ(DMPlexDistributeField(dm, sf, section, *v, dist_section, dist_v)); 38 for (p = 0; p < size; ++p) { 39 if (p == rank) { 40 CHKERRQ(PetscObjectViewFromOptions((PetscObject) dist_section, NULL, "-rd2_section_view")); 41 CHKERRQ(PetscObjectViewFromOptions((PetscObject) dist_v, NULL, "-rd2_vec_view")); 42 } 43 CHKERRQ(PetscBarrier((PetscObject) dm)); 44 CHKERRQ(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 45 } 46 47 CHKERRQ(PetscSectionDestroy(&dist_section)); 48 CHKERRQ(DMDestroy(&dist_v_dm)); 49 50 CHKERRQ(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 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 68 69 /* cells */ 70 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &start_cell, &end_cell)); 71 CHKERRQ(VecGetDM(cell_geom, &cell_dm)); 72 CHKERRQ(DMGetLocalSection(cell_dm, &cell_section)); 73 CHKERRQ(VecGetArrayRead(cell_geom, &cell_array)); 74 75 for (c = start_cell; c < end_cell; ++c) { 76 const PetscFVCellGeom *geom; 77 CHKERRQ(PetscSectionGetOffset(cell_section, c, &offset)); 78 geom = (PetscFVCellGeom*)&cell_array[offset]; 79 CHKERRQ(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "rank %d c %D 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 CHKERRQ(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL)); 82 CHKERRQ(VecRestoreArrayRead(cell_geom, &cell_array)); 83 84 /* faces */ 85 CHKERRQ(DMPlexGetHeightStratum(dm, 1, &start_face, &end_face)); 86 CHKERRQ(VecGetDM(face_geom, &face_dm)); 87 CHKERRQ(DMGetLocalSection(face_dm, &face_section)); 88 CHKERRQ(VecGetArrayRead(face_geom, &face_array)); 89 for (f = start_face; f < end_face; ++f) { 90 CHKERRQ(DMPlexGetSupport(dm, f, &cells)); 91 CHKERRQ(DMPlexGetSupportSize(dm, f, &supportSize)); 92 if (supportSize > 1) { 93 CHKERRQ(PetscSectionGetOffset(face_section, f, &offset)); 94 CHKERRQ(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "rank %d f %D cells %D,%D 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 CHKERRQ(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL)); 98 CHKERRQ(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 PetscErrorCode ierr; 110 111 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 112 CHKERRQ(DMCreate(PETSC_COMM_WORLD, &dm)); 113 CHKERRQ(DMSetType(dm, DMPLEX)); 114 CHKERRQ(DMSetFromOptions(dm)); 115 CHKERRQ(DMViewFromOptions(dm, NULL, "-dm_view")); 116 117 CHKERRQ(DMPlexComputeGeometryFVM(dm, &cell_geom, &face_geom)); 118 CHKERRQ(dm_view_geometry(dm, cell_geom, face_geom)); 119 120 /* redistribute */ 121 CHKERRQ(DMPlexGetPartitioner(dm, &part)); 122 CHKERRQ(PetscPartitionerSetFromOptions(part)); 123 CHKERRQ(PetscOptionsGetInt(NULL, NULL, "-overlap2", &overlap2, NULL)); 124 CHKERRQ(DMPlexDistribute(dm, overlap2, &redist_sf, &redist_dm)); 125 if (redist_dm) { 126 CHKERRQ(redistribute_vec(redist_dm, redist_sf, &cell_geom)); 127 CHKERRQ(redistribute_vec(redist_dm, redist_sf, &face_geom)); 128 CHKERRQ(PetscObjectViewFromOptions((PetscObject) redist_sf, NULL, "-rd2_sf_view")); 129 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "redistributed:\n")); 130 CHKERRQ(dm_view_geometry(redist_dm, cell_geom, face_geom)); 131 } 132 133 CHKERRQ(VecDestroy(&cell_geom)); 134 CHKERRQ(VecDestroy(&face_geom)); 135 CHKERRQ(PetscSFDestroy(&redist_sf)); 136 CHKERRQ(DMDestroy(&redist_dm)); 137 CHKERRQ(DMDestroy(&dm)); 138 ierr = PetscFinalize(); 139 return ierr; 140 } 141 142 /*TEST 143 144 test: 145 suffix: 0 146 nsize: 3 147 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 148 149 TEST*/ 150