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