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;
  8:   Vec          u;
  9:   PetscSection section;
 10:   PetscViewer  viewer;
 11:   PetscInt     dim, numFields, numBC, i;
 12:   PetscInt     numComp[3];
 13:   PetscInt     numDof[12];
 14:   PetscInt     bcField[1];
 15:   IS           bcPointIS[1];

 17:   PetscFunctionBeginUser;
 18:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 19:   /* Create a mesh */
 20:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 21:   PetscCall(DMSetType(dm, DMPLEX));
 22:   PetscCall(DMSetFromOptions(dm));
 23:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
 24:   PetscCall(DMGetDimension(dm, &dim));
 25:   /* Create a scalar field u, a vector field v, and a surface vector field w */
 26:   numFields  = 3;
 27:   numComp[0] = 1;
 28:   numComp[1] = dim;
 29:   numComp[2] = dim - 1;
 30:   for (i = 0; i < numFields * (dim + 1); ++i) numDof[i] = 0;
 31:   /* Let u be defined on vertices */
 32:   numDof[0 * (dim + 1) + 0] = 1;
 33:   /* Let v be defined on cells */
 34:   numDof[1 * (dim + 1) + dim] = dim;
 35:   /* Let w be defined on faces */
 36:   numDof[2 * (dim + 1) + dim - 1] = dim - 1;
 37:   /* Setup boundary conditions */
 38:   numBC = 1;
 39:   /* Prescribe a Dirichlet condition on u on the boundary
 40:        Label "marker" is made by the mesh creation routine */
 41:   bcField[0] = 0;
 42:   PetscCall(DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]));
 43:   /* Create a PetscSection with this data layout */
 44:   PetscCall(DMSetNumFields(dm, numFields));
 45:   PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section));
 46:   PetscCall(ISDestroy(&bcPointIS[0]));
 47:   /* Name the Field variables */
 48:   PetscCall(PetscSectionSetFieldName(section, 0, "u"));
 49:   PetscCall(PetscSectionSetFieldName(section, 1, "v"));
 50:   PetscCall(PetscSectionSetFieldName(section, 2, "w"));
 51:   PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD));
 52:   /* Tell the DM to use this data layout */
 53:   PetscCall(DMSetLocalSection(dm, section));
 54:   /* Create a Vec with this layout and view it */
 55:   PetscCall(DMGetGlobalVector(dm, &u));
 56:   PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
 57:   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
 58:   PetscCall(PetscViewerFileSetName(viewer, "sol.vtu"));
 59:   PetscCall(VecView(u, viewer));
 60:   PetscCall(PetscViewerDestroy(&viewer));
 61:   PetscCall(DMRestoreGlobalVector(dm, &u));
 62:   /* Cleanup */
 63:   PetscCall(PetscSectionDestroy(&section));
 64:   PetscCall(DMDestroy(&dm));
 65:   PetscCall(PetscFinalize());
 66:   return 0;
 67: }

 69: /*TEST

 71:   test:
 72:     suffix: 0
 73:     requires: triangle
 74:     args: -info :~sys,mat
 75:   test:
 76:     suffix: 1
 77:     requires: ctetgen
 78:     args: -dm_plex_dim 3 -info :~sys,mat

 80: TEST*/