Hi, I'm trying to get DMPlex to solve and display with Dirichlet boundary conditions. My global and local solution vectors are of different sizes. The length difference is the number of nodes read from the nodeset read in the exodus file. I am able to project the initial solution/boundary data that I want with with VecSetValuesSection. The initial condition looks ok with the global solution in the viewer. However, I seem to be losing the dirichlet boundary data on the solve ? I'm trying to use the suggested global to local scattering on the output. { Vec localX; PetscSection section; PetscInt vStart, vEnd, v ; DMGetLocalVector(dm, &localX); DMGetDefaultSection(dm, §ion); DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); for (v = vStart; v < vEnd; ++v) { PetscScalar bcData = 101.0; VecSetValuesSection(localX, section, v, &bcData , INSERT_BC_VALUES); } DMRestoreLocalVector(dm, &localX); DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, u); DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, u); PetscViewer viewer; char vtkfilename[PETSC_MAX_PATH_LEN] = "initialcondition.vtk"; ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,vtkfilename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); ierr = VecView(u,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr) } . . ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); . . if (write_output) { PetscViewer viewer; const char *name; Vec lv; DMGetLocalVector(dm, &lv); ierr = PetscObjectGetName((PetscObject) u,&name);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) lv, name);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr); char vtkfilename[PETSC_MAX_PATH_LEN] = "solution.vtk"; ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,vtkfilename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); ierr = VecView(lv,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); ierr = DMRestoreLocalVector(dm, &lv);CHKERRQ(ierr); } On Wed, Jan 29, 2014 at 2:04 PM, Geoffrey Irving <[email protected]> wrote:
On Wed, Jan 29, 2014 at 12:00 PM, Jed Brown <[email protected]> wrote:
Local vectors are supposed to be just that: local. VecViewing a local vector and expecting it to be parallel is perverse. So we need a real interface.
Placing 1.0 on the diagonal (and don't assemble into those rows and columns) is the common way to deal with Dirichlet boundary nodes. See ex48 for one example. I have written about this in a few places; I can find the more complete description when I have a keyboard.
I'll look forward to the improved interface. For better or worse, I'd like to be able to view simulations in the near future, so it looks like I have to go with Matt's perverse but working version now.
Geoffrey
On January 29, 2014 12:53:26 PM MST, Geoffrey Irving <[email protected]> wrote:
On Wed, Jan 29, 2014 at 11:36 AM, Jed Brown <[email protected]> wrote:
Matt's sample code doesn't set it either. We need to fix (and the *only* acceptable fix if that VecView always does the right thing, because we have to be able to call it in analysis settings that know nothing about your discretization).
Matt's sample code doesn't set it either, but for Matt's sample code I know where to insert the one line call to DMPlexProjectFunctionLocal. For your version I never have explicit access to the local vector, so I can't insert the fix.
The problem is that some vectors reside in a homogeneous space (e.g. increments and eigenvectors) while others reside in the inhomogeneous space (solutions). We can add a flag or BC attribute on the vector to this effect, but this (and slip conditions) was the issue that led me to conclude that removing boundary nodes was mostly a false economy.
To leave the boundary conditions in, we would need efficient support for a very large, very sparse MatNullSpace. This is doable via shells, but is it easy to do in a way that doesn't interfere with the user's other null spaces?
Geoffrey