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
redistribute_vec(DM dist_dm,PetscSF sf,Vec * v)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) PetscCall(PetscObjectViewFromOptions((PetscObject)*v, NULL, "-rd_vec_view"));
33 PetscCall(PetscBarrier((PetscObject)dm));
34 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
35 }
36 PetscCall(DMPlexDistributeField(dm, sf, section, *v, dist_section, dist_v));
37 for (p = 0; p < size; ++p) {
38 if (p == rank) {
39 PetscCall(PetscObjectViewFromOptions((PetscObject)dist_section, NULL, "-rd2_section_view"));
40 PetscCall(PetscObjectViewFromOptions((PetscObject)dist_v, NULL, "-rd2_vec_view"));
41 }
42 PetscCall(PetscBarrier((PetscObject)dm));
43 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
44 }
45
46 PetscCall(PetscSectionDestroy(&dist_section));
47 PetscCall(DMDestroy(&dist_v_dm));
48
49 PetscCall(VecDestroy(v));
50 *v = dist_v;
51 PetscFunctionReturn(PETSC_SUCCESS);
52 }
53
dm_view_geometry(DM dm,Vec cell_geom,Vec face_geom)54 static PetscErrorCode dm_view_geometry(DM dm, Vec cell_geom, Vec face_geom)
55 {
56 DM cell_dm, face_dm;
57 PetscSection cell_section, face_section;
58 const PetscScalar *cell_array, *face_array;
59 const PetscInt *cells;
60 PetscInt c, start_cell, end_cell;
61 PetscInt f, start_face, end_face;
62 PetscInt supportSize, offset;
63 PetscMPIInt rank;
64
65 PetscFunctionBegin;
66 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
67
68 /* cells */
69 PetscCall(DMPlexGetHeightStratum(dm, 0, &start_cell, &end_cell));
70 PetscCall(VecGetDM(cell_geom, &cell_dm));
71 PetscCall(DMGetLocalSection(cell_dm, &cell_section));
72 PetscCall(VecGetArrayRead(cell_geom, &cell_array));
73
74 for (c = start_cell; c < end_cell; ++c) {
75 const PetscFVCellGeom *geom;
76 PetscCall(PetscSectionGetOffset(cell_section, c, &offset));
77 geom = (PetscFVCellGeom *)&cell_array[offset];
78 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));
79 }
80 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
81 PetscCall(VecRestoreArrayRead(cell_geom, &cell_array));
82
83 /* faces */
84 PetscCall(DMPlexGetHeightStratum(dm, 1, &start_face, &end_face));
85 PetscCall(VecGetDM(face_geom, &face_dm));
86 PetscCall(DMGetLocalSection(face_dm, &face_section));
87 PetscCall(VecGetArrayRead(face_geom, &face_array));
88 for (f = start_face; f < end_face; ++f) {
89 PetscCall(DMPlexGetSupport(dm, f, &cells));
90 PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
91 if (supportSize > 1) {
92 PetscCall(PetscSectionGetOffset(face_section, f, &offset));
93 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])));
94 }
95 }
96 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
97 PetscCall(VecRestoreArrayRead(cell_geom, &cell_array));
98 PetscFunctionReturn(PETSC_SUCCESS);
99 }
100
main(int argc,char ** argv)101 int main(int argc, char **argv)
102 {
103 DM dm, redist_dm;
104 PetscPartitioner part;
105 PetscSF redist_sf;
106 Vec cell_geom, face_geom;
107 PetscInt overlap2 = 1;
108
109 PetscFunctionBeginUser;
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