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