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