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