Actual source code: ex1.c

      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:   PetscInt     dim, numFields, numBC;
     11:   PetscInt     numComp[3];
     12:   PetscInt     numDof[12];
     13:   PetscInt     bcField[1];
     14:   IS           bcPointIS[1];

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

     66: /*TEST

     68:   test:
     69:     suffix: 0
     70:     requires: triangle
     71:     args: -mysection_view -vec_view vtk:sol.vtu -info :~sys,mat

     73:   test:
     74:     suffix: 1
     75:     requires: ctetgen
     76:     args: -dm_plex_dim 3 -mysection_view -info :~sys,mat

     78:   test:
     79:     suffix: filter_0
     80:     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -dm_plex_label_lower 0,1,2,3,4,6 \
     81:       -dm_refine 1 -dm_plex_transform_type transform_filter -dm_plex_transform_active lower \
     82:       -dm_view

     84:   test:
     85:     suffix: submesh_0
     86:     requires: triangle
     87:     args: -dm_plex_label_middle 9,12,15,23,31 -dm_plex_submesh middle -dm_view

     89: TEST*/