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 CHKERRQ(PetscInitialize(&argc, &argv, NULL,help)); 18 /* Create a mesh */ 19 CHKERRQ(DMCreate(PETSC_COMM_WORLD, &dm)); 20 CHKERRQ(DMSetType(dm, DMPLEX)); 21 CHKERRQ(DMSetFromOptions(dm)); 22 CHKERRQ(DMViewFromOptions(dm, NULL, "-dm_view")); 23 CHKERRQ(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 (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 CHKERRQ(DMGetStratumIS(dm, "marker", 1, &bcPointIS[0])); 42 /* Create a PetscSection with this data layout */ 43 CHKERRQ(DMSetNumFields(dm, numFields)); 44 CHKERRQ(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion)); 45 CHKERRQ(ISDestroy(&bcPointIS[0])); 46 /* Name the Field variables */ 47 CHKERRQ(PetscSectionSetFieldName(section, 0, "u")); 48 CHKERRQ(PetscSectionSetFieldName(section, 1, "v")); 49 CHKERRQ(PetscSectionSetFieldName(section, 2, "w")); 50 CHKERRQ(PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD)); 51 /* Tell the DM to use this data layout */ 52 CHKERRQ(DMSetLocalSection(dm, section)); 53 /* Create a Vec with this layout and view it */ 54 CHKERRQ(DMGetGlobalVector(dm, &u)); 55 CHKERRQ(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 56 CHKERRQ(PetscViewerSetType(viewer, PETSCVIEWERVTK)); 57 CHKERRQ(PetscViewerFileSetName(viewer, "sol.vtu")); 58 CHKERRQ(VecView(u, viewer)); 59 CHKERRQ(PetscViewerDestroy(&viewer)); 60 CHKERRQ(DMRestoreGlobalVector(dm, &u)); 61 /* Cleanup */ 62 CHKERRQ(PetscSectionDestroy(§ion)); 63 CHKERRQ(DMDestroy(&dm)); 64 CHKERRQ(PetscFinalize()); 65 return 0; 66 } 67 68 /*TEST 69 70 test: 71 suffix: 0 72 requires: triangle 73 args: -info :~sys,mat 74 test: 75 suffix: 1 76 requires: ctetgen 77 args: -dm_plex_dim 3 -info :~sys,mat 78 79 TEST*/ 80