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 PetscInt dim, numFields, numBC; 11 PetscInt numComp[3]; 12 PetscInt numDof[12]; 13 PetscInt bcField[1]; 14 IS bcPointIS[1]; 15 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, §ion)); 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(§ion)); 61 PetscCall(DMDestroy(&dm)); 62 PetscCall(PetscFinalize()); 63 return 0; 64 } 65 66 /*TEST 67 68 test: 69 suffix: 0 70 requires: triangle 71 args: -mysection_view -vec_view vtk:sol.vtu -info :~sys,mat 72 73 test: 74 suffix: 1 75 requires: ctetgen 76 args: -dm_plex_dim 3 -mysection_view -info :~sys,mat 77 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 83 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 88 89 TEST*/ 90