xref: /petsc/src/dm/impls/plex/tutorials/ex1.c (revision b122ec5aa1bd4469eb4e0673542fb7de3f411254)
1c4762a1bSJed Brown static char help[] = "Define a simple field over the mesh\n\n";
2c4762a1bSJed Brown 
3c4762a1bSJed Brown #include <petscdmplex.h>
4c4762a1bSJed Brown 
5c4762a1bSJed Brown int main(int argc, char **argv)
6c4762a1bSJed Brown {
730602db0SMatthew G. Knepley   DM             dm;
8c4762a1bSJed Brown   Vec            u;
9c4762a1bSJed Brown   PetscSection   section;
10c4762a1bSJed Brown   PetscViewer    viewer;
1130602db0SMatthew G. Knepley   PetscInt       dim, numFields, numBC, i;
12c4762a1bSJed Brown   PetscInt       numComp[3];
13c4762a1bSJed Brown   PetscInt       numDof[12];
14c4762a1bSJed Brown   PetscInt       bcField[1];
15c4762a1bSJed Brown   IS             bcPointIS[1];
16c4762a1bSJed Brown 
17*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscInitialize(&argc, &argv, NULL,help));
18c4762a1bSJed Brown   /* Create a mesh */
195f80ce2aSJacob Faibussowitsch   CHKERRQ(DMCreate(PETSC_COMM_WORLD, &dm));
205f80ce2aSJacob Faibussowitsch   CHKERRQ(DMSetType(dm, DMPLEX));
215f80ce2aSJacob Faibussowitsch   CHKERRQ(DMSetFromOptions(dm));
225f80ce2aSJacob Faibussowitsch   CHKERRQ(DMViewFromOptions(dm, NULL, "-dm_view"));
235f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetDimension(dm, &dim));
24c4762a1bSJed Brown   /* Create a scalar field u, a vector field v, and a surface vector field w */
25c4762a1bSJed Brown   numFields  = 3;
26c4762a1bSJed Brown   numComp[0] = 1;
27c4762a1bSJed Brown   numComp[1] = dim;
28c4762a1bSJed Brown   numComp[2] = dim-1;
29c4762a1bSJed Brown   for (i = 0; i < numFields*(dim+1); ++i) numDof[i] = 0;
30c4762a1bSJed Brown   /* Let u be defined on vertices */
31c4762a1bSJed Brown   numDof[0*(dim+1)+0]     = 1;
32c4762a1bSJed Brown   /* Let v be defined on cells */
33c4762a1bSJed Brown   numDof[1*(dim+1)+dim]   = dim;
34c4762a1bSJed Brown   /* Let w be defined on faces */
35c4762a1bSJed Brown   numDof[2*(dim+1)+dim-1] = dim-1;
36c4762a1bSJed Brown   /* Setup boundary conditions */
37c4762a1bSJed Brown   numBC = 1;
38c4762a1bSJed Brown   /* Prescribe a Dirichlet condition on u on the boundary
39c4762a1bSJed Brown        Label "marker" is made by the mesh creation routine */
40c4762a1bSJed Brown   bcField[0] = 0;
415f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]));
42c4762a1bSJed Brown   /* Create a PetscSection with this data layout */
435f80ce2aSJacob Faibussowitsch   CHKERRQ(DMSetNumFields(dm, numFields));
445f80ce2aSJacob Faibussowitsch   CHKERRQ(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section));
455f80ce2aSJacob Faibussowitsch   CHKERRQ(ISDestroy(&bcPointIS[0]));
46c4762a1bSJed Brown   /* Name the Field variables */
475f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionSetFieldName(section, 0, "u"));
485f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionSetFieldName(section, 1, "v"));
495f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionSetFieldName(section, 2, "w"));
505f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD));
51c4762a1bSJed Brown   /* Tell the DM to use this data layout */
525f80ce2aSJacob Faibussowitsch   CHKERRQ(DMSetLocalSection(dm, section));
53c4762a1bSJed Brown   /* Create a Vec with this layout and view it */
545f80ce2aSJacob Faibussowitsch   CHKERRQ(DMGetGlobalVector(dm, &u));
555f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
565f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerSetType(viewer, PETSCVIEWERVTK));
575f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerFileSetName(viewer, "sol.vtu"));
585f80ce2aSJacob Faibussowitsch   CHKERRQ(VecView(u, viewer));
595f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscViewerDestroy(&viewer));
605f80ce2aSJacob Faibussowitsch   CHKERRQ(DMRestoreGlobalVector(dm, &u));
61c4762a1bSJed Brown   /* Cleanup */
625f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscSectionDestroy(&section));
635f80ce2aSJacob Faibussowitsch   CHKERRQ(DMDestroy(&dm));
64*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscFinalize());
65*b122ec5aSJacob Faibussowitsch   return 0;
66c4762a1bSJed Brown }
67c4762a1bSJed Brown 
68c4762a1bSJed Brown /*TEST
69c4762a1bSJed Brown 
70c4762a1bSJed Brown   test:
71c4762a1bSJed Brown     suffix: 0
72c4762a1bSJed Brown     requires: triangle
7330602db0SMatthew G. Knepley     args: -info :~sys,mat
74c4762a1bSJed Brown   test:
75c4762a1bSJed Brown     suffix: 1
76c4762a1bSJed Brown     requires: ctetgen
7730602db0SMatthew G. Knepley     args: -dm_plex_dim 3 -info :~sys,mat
78c4762a1bSJed Brown 
79c4762a1bSJed Brown TEST*/
80