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);CHKERRQ(ierr); 31 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(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);CHKERRQ(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, dist_dm, redist_dm; 107 PetscPartitioner part; 108 PetscSF dist_sf, redist_sf; 109 Vec cell_geom, face_geom; 110 PetscInt overlap = 1, overlap2 = 1; 111 PetscMPIInt rank; 112 const char *filename = "gminc_1d.exo"; 113 PetscErrorCode ierr; 114 115 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 116 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr); 117 if (0) { 118 ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, PETSC_TRUE, &dm);CHKERRQ(ierr); 119 } else { 120 ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 3, PETSC_FALSE, NULL, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr); 121 } 122 ierr = DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);CHKERRQ(ierr); 123 124 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 125 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 126 ierr = PetscOptionsGetInt(NULL, NULL, "-overlap", &overlap, NULL);CHKERRQ(ierr); 127 ierr = DMPlexDistribute(dm, overlap, &dist_sf, &dist_dm);CHKERRQ(ierr); 128 if (dist_dm) { 129 ierr = DMDestroy(&dm);CHKERRQ(ierr); 130 dm = dist_dm; 131 } 132 133 ierr = DMPlexComputeGeometryFVM(dm, &cell_geom, &face_geom);CHKERRQ(ierr); 134 135 ierr = dm_view_geometry(dm, cell_geom, face_geom);CHKERRQ(ierr); 136 137 /* redistribute */ 138 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 139 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 140 ierr = PetscOptionsGetInt(NULL, NULL, "-overlap2", &overlap2, NULL);CHKERRQ(ierr); 141 ierr = DMPlexDistribute(dm, overlap2, &redist_sf, &redist_dm);CHKERRQ(ierr); 142 if (redist_dm) { 143 ierr = redistribute_vec(redist_dm, redist_sf, &cell_geom);CHKERRQ(ierr); 144 ierr = redistribute_vec(redist_dm, redist_sf, &face_geom);CHKERRQ(ierr); 145 ierr = PetscObjectViewFromOptions((PetscObject) redist_sf, NULL, "-rd2_sf_view");CHKERRQ(ierr); 146 ierr = PetscPrintf(PETSC_COMM_WORLD, "redistributed:\n");CHKERRQ(ierr); 147 ierr = dm_view_geometry(redist_dm, cell_geom, face_geom);CHKERRQ(ierr); 148 } 149 150 ierr = VecDestroy(&cell_geom);CHKERRQ(ierr); 151 ierr = VecDestroy(&face_geom);CHKERRQ(ierr); 152 ierr = PetscSFDestroy(&dist_sf);CHKERRQ(ierr); 153 ierr = PetscSFDestroy(&redist_sf);CHKERRQ(ierr); 154 ierr = DMDestroy(&redist_dm);CHKERRQ(ierr); 155 ierr = DMDestroy(&dm);CHKERRQ(ierr); 156 ierr = PetscFinalize(); 157 return ierr; 158 } 159 160 /*TEST 161 162 test: 163 suffix: 0 164 nsize: 3 165 args: -overlap 1 -overlap2 1 -dm_plex_box_faces 8,1,1 -petscpartitioner_type simple 166 167 TEST*/ 168