xref: /petsc/src/dm/impls/plex/tutorials/ex1.c (revision 327415f76d85372a4417cf1aaa14db707d4d6c04)
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*327415f7SBarry Smith   PetscFunctionBeginUser;
189566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, NULL,help));
19c4762a1bSJed Brown   /* Create a mesh */
209566063dSJacob Faibussowitsch   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
219566063dSJacob Faibussowitsch   PetscCall(DMSetType(dm, DMPLEX));
229566063dSJacob Faibussowitsch   PetscCall(DMSetFromOptions(dm));
239566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
249566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
25c4762a1bSJed Brown   /* Create a scalar field u, a vector field v, and a surface vector field w */
26c4762a1bSJed Brown   numFields  = 3;
27c4762a1bSJed Brown   numComp[0] = 1;
28c4762a1bSJed Brown   numComp[1] = dim;
29c4762a1bSJed Brown   numComp[2] = dim-1;
30c4762a1bSJed Brown   for (i = 0; i < numFields*(dim+1); ++i) numDof[i] = 0;
31c4762a1bSJed Brown   /* Let u be defined on vertices */
32c4762a1bSJed Brown   numDof[0*(dim+1)+0]     = 1;
33c4762a1bSJed Brown   /* Let v be defined on cells */
34c4762a1bSJed Brown   numDof[1*(dim+1)+dim]   = dim;
35c4762a1bSJed Brown   /* Let w be defined on faces */
36c4762a1bSJed Brown   numDof[2*(dim+1)+dim-1] = dim-1;
37c4762a1bSJed Brown   /* Setup boundary conditions */
38c4762a1bSJed Brown   numBC = 1;
39c4762a1bSJed Brown   /* Prescribe a Dirichlet condition on u on the boundary
40c4762a1bSJed Brown        Label "marker" is made by the mesh creation routine */
41c4762a1bSJed Brown   bcField[0] = 0;
429566063dSJacob Faibussowitsch   PetscCall(DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]));
43c4762a1bSJed Brown   /* Create a PetscSection with this data layout */
449566063dSJacob Faibussowitsch   PetscCall(DMSetNumFields(dm, numFields));
459566063dSJacob Faibussowitsch   PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section));
469566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&bcPointIS[0]));
47c4762a1bSJed Brown   /* Name the Field variables */
489566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetFieldName(section, 0, "u"));
499566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetFieldName(section, 1, "v"));
509566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetFieldName(section, 2, "w"));
519566063dSJacob Faibussowitsch   PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD));
52c4762a1bSJed Brown   /* Tell the DM to use this data layout */
539566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(dm, section));
54c4762a1bSJed Brown   /* Create a Vec with this layout and view it */
559566063dSJacob Faibussowitsch   PetscCall(DMGetGlobalVector(dm, &u));
569566063dSJacob Faibussowitsch   PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
579566063dSJacob Faibussowitsch   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
589566063dSJacob Faibussowitsch   PetscCall(PetscViewerFileSetName(viewer, "sol.vtu"));
599566063dSJacob Faibussowitsch   PetscCall(VecView(u, viewer));
609566063dSJacob Faibussowitsch   PetscCall(PetscViewerDestroy(&viewer));
619566063dSJacob Faibussowitsch   PetscCall(DMRestoreGlobalVector(dm, &u));
62c4762a1bSJed Brown   /* Cleanup */
639566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&section));
649566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dm));
659566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
66b122ec5aSJacob Faibussowitsch   return 0;
67c4762a1bSJed Brown }
68c4762a1bSJed Brown 
69c4762a1bSJed Brown /*TEST
70c4762a1bSJed Brown 
71c4762a1bSJed Brown   test:
72c4762a1bSJed Brown     suffix: 0
73c4762a1bSJed Brown     requires: triangle
7430602db0SMatthew G. Knepley     args: -info :~sys,mat
75c4762a1bSJed Brown   test:
76c4762a1bSJed Brown     suffix: 1
77c4762a1bSJed Brown     requires: ctetgen
7830602db0SMatthew G. Knepley     args: -dm_plex_dim 3 -info :~sys,mat
79c4762a1bSJed Brown 
80c4762a1bSJed Brown TEST*/
81