Actual source code: ex66.c

  1: static const char help[] = "Test for non-manifold interpolation";

  3: #include <petscdmplex.h>

  5: /*
  6:        3-------------7
  7:       /|            /|
  8:      / |           / |
  9:     /  |          /  |
 10:    1-------------5   |
 11:    |   |         |   |
 12:    |   |         |   |
 13:    |   |         |   |
 14:    |   |         |   |
 15:    z   4---------|---8
 16:    ^  /          |  /
 17:    | y           | /
 18:    |/            |/
 19:    2--->-x-------6-------------9
 20: */
 21: int main(int argc, char **argv)
 22: {
 23:   DM        dm, idm;
 24:   DMLabel   ctLabel;
 25:   PetscBool has_vtk = PETSC_FALSE;

 27:   // 9 vertices
 28:   // 1 edge
 29:   // 0 faces
 30:   // 1 volume
 31:   PetscInt num_points[4] = {9, 1, 0, 1};

 33:   // point 0 = hexahedron (defined by 8 vertices)
 34:   // points 1-9 = vertices
 35:   // point 10 = edged (defined by 2 vertices)
 36:   PetscInt cone_size[11] = {8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2};

 38:   // hexahedron defined by points
 39:   PetscInt    cones[11]             = {3, 4, 2, 1, 7, 5, 6, 8, 6, 9};
 40:   PetscInt    cone_orientations[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 41:   PetscScalar vertex_coords[3 * 9]  = {0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 2, 0, 0};

 43:   PetscFunctionBeginUser;
 44:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 45:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Output VTK?", "ex66");
 46:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-vtk", &has_vtk, NULL));
 47:   PetscOptionsEnd();

 49:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 50:   PetscCall(PetscObjectSetName((PetscObject)dm, "cubeline-fromdag"));
 51:   PetscCall(DMSetType(dm, DMPLEX));
 52:   PetscCall(DMSetDimension(dm, 3));

 54:   PetscCall(DMPlexCreateFromDAG(dm, 3, num_points, cone_size, cones, cone_orientations, vertex_coords));
 55:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

 57:   // TODO: make it work with a DM made from a msh file
 58:   // PetscCall(DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, "cube-line.msh", PETSC_FALSE, &dm));
 59:   // PetscCall(PetscObjectSetName((PetscObject)dm, "msh"));

 61:   // Must set cell types
 62:   PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
 63:   PetscCall(DMLabelSetValue(ctLabel, 0, DM_POLYTOPE_HEXAHEDRON));
 64:   PetscCall(DMLabelSetValue(ctLabel, 1, DM_POLYTOPE_POINT));
 65:   PetscCall(DMLabelSetValue(ctLabel, 2, DM_POLYTOPE_POINT));
 66:   PetscCall(DMLabelSetValue(ctLabel, 3, DM_POLYTOPE_POINT));
 67:   PetscCall(DMLabelSetValue(ctLabel, 4, DM_POLYTOPE_POINT));
 68:   PetscCall(DMLabelSetValue(ctLabel, 5, DM_POLYTOPE_POINT));
 69:   PetscCall(DMLabelSetValue(ctLabel, 6, DM_POLYTOPE_POINT));
 70:   PetscCall(DMLabelSetValue(ctLabel, 7, DM_POLYTOPE_POINT));
 71:   PetscCall(DMLabelSetValue(ctLabel, 8, DM_POLYTOPE_POINT));
 72:   PetscCall(DMLabelSetValue(ctLabel, 9, DM_POLYTOPE_POINT));
 73:   PetscCall(DMLabelSetValue(ctLabel, 10, DM_POLYTOPE_SEGMENT));

 75:   // interpolate (make sure to use -interp_dm_plex_stratify_celltype)
 76:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, "interp_"));
 77:   PetscCall(DMPlexInterpolate(dm, &idm));
 78:   PetscCall(DMDestroy(&dm));
 79:   dm = idm;

 81:   PetscCall(DMSetFromOptions(dm));
 82:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

 84:   if (has_vtk) {
 85:     PetscViewer viewer;
 86:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
 87:     PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
 88:     PetscCall(PetscViewerFileSetName(viewer, "ex66.vtk"));
 89:     PetscCall(DMView(dm, viewer));
 90:     PetscCall(PetscViewerDestroy(&viewer));
 91:   }

 93:   PetscCall(DMDestroy(&dm));
 94:   PetscCall(PetscFinalize());
 95:   return 0;
 96: }

 98: /*TEST
 99:   test:
100:     suffix: 0
101:     args: -interp_dm_plex_stratify_celltype -dm_view ::ascii_info_detail -interp_dm_view ::ascii_info_detail

103: TEST*/