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, dmDist = NULL; 8 Vec u; 9 PetscSection section; 10 PetscViewer viewer; 11 PetscInt dim = 2, numFields, numBC, i; 12 PetscInt numComp[3]; 13 PetscInt numDof[12]; 14 PetscInt bcField[1]; 15 IS bcPointIS[1]; 16 PetscBool interpolate = PETSC_TRUE; 17 PetscErrorCode ierr; 18 19 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 20 ierr = PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);CHKERRQ(ierr); 21 /* Create a mesh */ 22 ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, NULL, NULL, NULL, NULL, interpolate, &dm);CHKERRQ(ierr); 23 /* Distribute mesh over processes */ 24 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);CHKERRQ(ierr); 25 if (dmDist) {ierr = DMDestroy(&dm);CHKERRQ(ierr); dm = dmDist;} 26 /* Create a scalar field u, a vector field v, and a surface vector field w */ 27 numFields = 3; 28 numComp[0] = 1; 29 numComp[1] = dim; 30 numComp[2] = dim-1; 31 for (i = 0; i < numFields*(dim+1); ++i) numDof[i] = 0; 32 /* Let u be defined on vertices */ 33 numDof[0*(dim+1)+0] = 1; 34 /* Let v be defined on cells */ 35 numDof[1*(dim+1)+dim] = dim; 36 /* Let w be defined on faces */ 37 numDof[2*(dim+1)+dim-1] = dim-1; 38 /* Setup boundary conditions */ 39 numBC = 1; 40 /* Prescribe a Dirichlet condition on u on the boundary 41 Label "marker" is made by the mesh creation routine */ 42 bcField[0] = 0; 43 ierr = DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]);CHKERRQ(ierr); 44 /* Create a PetscSection with this data layout */ 45 ierr = DMSetNumFields(dm, numFields);CHKERRQ(ierr); 46 ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion);CHKERRQ(ierr); 47 ierr = ISDestroy(&bcPointIS[0]);CHKERRQ(ierr); 48 /* Name the Field variables */ 49 ierr = PetscSectionSetFieldName(section, 0, "u");CHKERRQ(ierr); 50 ierr = PetscSectionSetFieldName(section, 1, "v");CHKERRQ(ierr); 51 ierr = PetscSectionSetFieldName(section, 2, "w");CHKERRQ(ierr); 52 ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 53 /* Tell the DM to use this data layout */ 54 ierr = DMSetLocalSection(dm, section);CHKERRQ(ierr); 55 /* Create a Vec with this layout and view it */ 56 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 57 ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr); 58 ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr); 59 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr); 60 ierr = PetscViewerFileSetName(viewer, "sol.vtk");CHKERRQ(ierr); 61 ierr = VecView(u, viewer);CHKERRQ(ierr); 62 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 63 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 64 /* Cleanup */ 65 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 66 ierr = DMDestroy(&dm);CHKERRQ(ierr); 67 ierr = PetscFinalize(); 68 return ierr; 69 } 70 71 /*TEST 72 73 test: 74 suffix: 0 75 requires: triangle 76 args: -info :~sys 77 test: 78 suffix: 1 79 requires: ctetgen 80 args: -dim 3 -info :~sys 81 82 TEST*/ 83