Actual source code: ex1.c
1: static char help[] = "Define a simple field over the mesh\n\n";
3: #include <petscdmplex.h>
5: int main(int argc, char **argv)
6: {
7: DM dm;
8: Vec u;
9: PetscSection section;
10: PetscViewer viewer;
11: PetscInt dim, numFields, numBC, i;
12: PetscInt numComp[3];
13: PetscInt numDof[12];
14: PetscInt bcField[1];
15: IS bcPointIS[1];
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
19: /* Create a mesh */
20: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
21: PetscCall(DMSetType(dm, DMPLEX));
22: PetscCall(DMSetFromOptions(dm));
23: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
24: PetscCall(DMGetDimension(dm, &dim));
25: /* Create a scalar field u, a vector field v, and a surface vector field w */
26: numFields = 3;
27: numComp[0] = 1;
28: numComp[1] = dim;
29: numComp[2] = dim - 1;
30: for (i = 0; i < numFields * (dim + 1); ++i) numDof[i] = 0;
31: /* Let u be defined on vertices */
32: numDof[0 * (dim + 1) + 0] = 1;
33: /* Let v be defined on cells */
34: numDof[1 * (dim + 1) + dim] = dim;
35: /* Let w be defined on faces */
36: numDof[2 * (dim + 1) + dim - 1] = dim - 1;
37: /* Setup boundary conditions */
38: numBC = 1;
39: /* Prescribe a Dirichlet condition on u on the boundary
40: Label "marker" is made by the mesh creation routine */
41: bcField[0] = 0;
42: PetscCall(DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]));
43: /* Create a PetscSection with this data layout */
44: PetscCall(DMSetNumFields(dm, numFields));
45: PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion));
46: PetscCall(ISDestroy(&bcPointIS[0]));
47: /* Name the Field variables */
48: PetscCall(PetscSectionSetFieldName(section, 0, "u"));
49: PetscCall(PetscSectionSetFieldName(section, 1, "v"));
50: PetscCall(PetscSectionSetFieldName(section, 2, "w"));
51: PetscCall(PetscSectionViewFromOptions(section, NULL, "-mysection_view"));
52: /* Tell the DM to use this data layout */
53: PetscCall(DMSetLocalSection(dm, section));
54: /* Create a Vec with this layout and view it */
55: PetscCall(DMGetGlobalVector(dm, &u));
56: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
57: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
58: PetscCall(PetscViewerFileSetName(viewer, "sol.vtu"));
59: PetscCall(VecView(u, viewer));
60: PetscCall(PetscViewerDestroy(&viewer));
61: PetscCall(DMRestoreGlobalVector(dm, &u));
62: /* Cleanup */
63: PetscCall(PetscSectionDestroy(§ion));
64: PetscCall(DMDestroy(&dm));
65: PetscCall(PetscFinalize());
66: return 0;
67: }
69: /*TEST
71: test:
72: suffix: 0
73: requires: triangle
74: args: -mysection_view -info :~sys,mat
76: test:
77: suffix: 1
78: requires: ctetgen
79: args: -dm_plex_dim 3 -mysection_view -info :~sys,mat
81: test:
82: suffix: filter_0
83: args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -dm_plex_label_lower 0,1,2,3,4,6 \
84: -dm_refine 1 -dm_plex_transform_type transform_filter -dm_plex_transform_active lower \
85: -dm_view
87: test:
88: suffix: submesh_0
89: requires: triangle
90: args: -dm_plex_label_middle 9,12,15,23,31 -dm_plex_submesh middle -dm_view
92: TEST*/