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