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, &section);
 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*/