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, PETSC_TRUE, NULL, NULL, NULL, NULL, 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 w 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:   DMSetNumFields(dm, numFields);
 46:   DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section);
 47:   ISDestroy(&bcPointIS[0]);
 48:   /* Name the Field variables */
 49:   PetscSectionSetFieldName(section, 0, "u");
 50:   PetscSectionSetFieldName(section, 1, "v");
 51:   PetscSectionSetFieldName(section, 2, "w");
 52:   PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);
 53:   /* Tell the DM to use this data layout */
 54:   DMSetLocalSection(dm, section);
 55:   /* Create a Vec with this layout and view it */
 56:   DMGetGlobalVector(dm, &u);
 57:   PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
 58:   PetscViewerSetType(viewer, PETSCVIEWERVTK);
 59:   PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK);
 60:   PetscViewerFileSetName(viewer, "sol.vtk");
 61:   VecView(u, viewer);
 62:   PetscViewerDestroy(&viewer);
 63:   DMRestoreGlobalVector(dm, &u);
 64:   /* Cleanup */
 65:   PetscSectionDestroy(&section);
 66:   DMDestroy(&dm);
 67:   PetscFinalize();
 68:   return ierr;
 69: }

 71: /*TEST

 73:   test:
 74:     suffix: 0
 75:     requires: triangle
 76:     args: -info :~sys
 77:   test:
 78:     suffix: 1
 79:     requires: ctetgen
 80:     args: -dim 3 -info :~sys

 82: TEST*/