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