1 static char help[] = "Define a simple field over the mesh\n\n"; 2 3 #include <petscdmplex.h> 4 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]; 16 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, §ion)); 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(PetscSectionViewFromOptions(section, NULL, "-mysection_view")); 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(§ion)); 64 PetscCall(DMDestroy(&dm)); 65 PetscCall(PetscFinalize()); 66 return 0; 67 } 68 69 /*TEST 70 71 test: 72 suffix: 0 73 requires: triangle 74 args: -mysection_view -info :~sys,mat 75 76 test: 77 suffix: 1 78 requires: ctetgen 79 args: -dm_plex_dim 3 -mysection_view -info :~sys,mat 80 81 test: 82 suffix: filter_0 83 args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -dm_plex_label_lower 0,1,2,3,4,6 \ 84 -dm_refine 1 -dm_plex_transform_type transform_filter -dm_plex_transform_active lower \ 85 -dm_view 86 87 test: 88 suffix: submesh_0 89 requires: triangle 90 args: -dm_plex_label_middle 9,12,15,23,31 -dm_plex_submesh middle -dm_view 91 92 TEST*/ 93