Actual source code: ex1.c
petsc-3.7.7 2017-09-25
1: static char help[] = "Define a simple field over the mesh\n\n";
3: #include <petscdmplex.h>
7: int main(int argc, char **argv)
8: {
9: DM dm, dmDist = NULL;
10: Vec u;
11: PetscSection section;
12: PetscViewer viewer;
13: PetscInt dim = 2, numFields, numBC, i;
14: PetscInt numComp[3];
15: PetscInt numDof[12];
16: PetscInt bcField[1];
17: IS bcPointIS[1];
18: PetscBool interpolate = PETSC_TRUE;
21: PetscInitialize(&argc, &argv, NULL, help);
22: PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);
23: /* Create a mesh */
24: DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, interpolate, &dm);
25: /* Distribute mesh over processes */
26: DMPlexDistribute(dm, 0, NULL, &dmDist);
27: if (dmDist) {DMDestroy(&dm); dm = dmDist;}
28: /* Create a scalar field u, a vector field v, and a surface vector field w */
29: numFields = 3;
30: numComp[0] = 1;
31: numComp[1] = dim;
32: numComp[2] = dim-1;
33: for (i = 0; i < numFields*(dim+1); ++i) numDof[i] = 0;
34: /* Let u be defined on vertices */
35: numDof[0*(dim+1)+0] = 1;
36: /* Let v be defined on cells */
37: numDof[1*(dim+1)+dim] = dim;
38: /* Let v be defined on faces */
39: numDof[2*(dim+1)+dim-1] = dim-1;
40: /* Setup boundary conditions */
41: numBC = 1;
42: /* Prescribe a Dirichlet condition on u on the boundary
43: Label "marker" is made by the mesh creation routine */
44: bcField[0] = 0;
45: DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]);
46: /* Create a PetscSection with this data layout */
47: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion);
48: ISDestroy(&bcPointIS[0]);
49: /* Name the Field variables */
50: PetscSectionSetFieldName(section, 0, "u");
51: PetscSectionSetFieldName(section, 1, "v");
52: PetscSectionSetFieldName(section, 2, "w");
53: PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);
54: /* Tell the DM to use this data layout */
55: DMSetDefaultSection(dm, section);
56: /* Create a Vec with this layout and view it */
57: DMGetGlobalVector(dm, &u);
58: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
59: PetscViewerSetType(viewer, PETSCVIEWERVTK);
60: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK);
61: PetscViewerFileSetName(viewer, "sol.vtk");
62: VecView(u, viewer);
63: PetscViewerDestroy(&viewer);
64: DMRestoreGlobalVector(dm, &u);
65: /* Cleanup */
66: PetscSectionDestroy(§ion);
67: DMDestroy(&dm);
68: PetscFinalize();
69: return 0;
70: }