In my code, I apply boundary conditions to the local vector before it is used for assembling the residual and jacobian. I use VecSetValuesSection with mode= INSERT_BC_VALUES
I also do the boundary condition's insert on the local vector before it is used for vecview. On Wed, Jan 29, 2014 at 12:00 PM, Geoffrey Irving <[email protected]> wrote: > On Tue, Jan 28, 2014 at 8:57 PM, Geoffrey Irving <[email protected]> wrote: >> On Tue, Jan 28, 2014 at 8:24 PM, Jed Brown <[email protected]> wrote: >>> Matthew Knepley <[email protected]> writes: >>> >>>> On Tue, Jan 28, 2014 at 9:41 PM, Jed Brown <[email protected]> wrote: >>>> >>>>> Matthew Knepley <[email protected]> writes: >>>>> > In ex62: >>>>> > >>>>> > if (user.runType == RUN_FULL) { >>>>> > PetscViewer viewer; >>>>> > Vec uLocal; >>>>> > const char *name; >>>>> > >>>>> > ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr); >>>>> > ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr); >>>>> > ierr = PetscViewerSetFormat(viewer, >>>>> > PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr); >>>>> > ierr = PetscViewerFileSetName(viewer, "ex62_sol.vtk");CHKERRQ(ierr); >>>>> >>>>> Please don't do this. Just set the type to PETSCVIEWERVTK, set the file >>>>> name to some *.vtu; the VTK viewer will infer output format from >>>>> extension and write the VTU XML format with binary-appended data in a >>>>> fairly fast and memory-scalable way. Then >>>>> >>>>> VecView(u,viewer); >>>> >>>> I think Jed is wrong here. You need to be using local vectors since >>>> they have BC and constraints in them. Please explain how your VTU >>>> format will magically put them in. >>> >>> We made this work when I threw a temper tantrum while we were writing TS >>> ex11.c (late 2012). >>> >>> PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer) >>> { >>> DM dm; >>> PetscBool isvtk; >>> PetscErrorCode ierr; >>> >>> PetscFunctionBegin; >>> ierr = VecGetDM(v, &dm);CHKERRQ(ierr); >>> if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, >>> "Vector not generated from a DM"); >>> ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, >>> &isvtk);CHKERRQ(ierr); >>> if (isvtk) { >>> Vec locv; >>> const char *name; >>> >>> ierr = DMGetLocalVector(dm, &locv);CHKERRQ(ierr); >>> ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr); >>> ierr = PetscObjectSetName((PetscObject) locv, name);CHKERRQ(ierr); >>> ierr = DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr); >>> ierr = DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr); >>> ierr = VecView_Plex_Local(locv, viewer);CHKERRQ(ierr); >>> ierr = DMRestoreLocalVector(dm, &locv);CHKERRQ(ierr); >> >> Which of those lines adds the correct boundary condition values? > > It looks like one of the following is true: > > 1. VecView on a DMPlex global vector does not insert correct boundary > conditions. > 2. Scalar, residual, and jacobian integration are all performing an > extra boundary condition enforcement step. > > I'm guessing it's (1). If so, does that mean I have to do it Matt's > more verbose way? > > Geoffrey
