<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Mon, Mar 31, 2014 at 9:27 AM, David Fuentes <span dir="ltr"><<a href="mailto:fuentesdt@gmail.com" target="_blank">fuentesdt@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Hi,</div><div><br></div><div>I'm trying to get DMPlex to solve and display with Dirichlet boundary conditions.</div>
<div><br></div><div>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.</div>
<div>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</div><div>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.</div>
</div></blockquote><div><br></div><div>I have just worked all this out in 'next', and most of it is in 'master' now. I will talk through the process</div><div><br></div><div>1) Suppose you view a global vector, either VTK or HDF5</div>
<div><br></div><div>2) PETSc maps this to a local vector, and calls the local version, which is the same as viewing a local vector</div><div><br></div><div>3) PETSc will insert the correct boundary values (via DMPlexInsertBoundaryValuesFEM) if you specify them</div>
<div> using the Section interface + DMAddBoundary(). Otherwise, you can view a local vector with the correct boundary values.</div><div><br></div><div>4) A scatter is created to a global vector including the boundary values</div>
<div><br></div><div>5) This expanded global vector is passes to the normal HDF5/VTK routines</div><div><br></div><div>It is still a little rough, but it should not be much effort to get it working the way you want.</div><div>
<br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div> { </div>
<div> Vec localX; </div>
<div>PetscSection section; </div>
<div> PetscInt vStart, vEnd, v ;</div><div><br></div><div> DMGetLocalVector(dm, &localX);</div><div> DMGetDefaultSection(dm, §ion);</div><div><br></div><div> DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);</div>
<div> for (v = vStart; v < vEnd; ++v) {</div><div> PetscScalar bcData = 101.0;</div><div> VecSetValuesSection(localX, section, v, &bcData , INSERT_BC_VALUES);</div><div> }</div><div> DMRestoreLocalVector(dm, &localX);</div>
<div> DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, u);</div><div> DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, u);</div><div><br></div><div> PetscViewer viewer;</div><div> char vtkfilename[PETSC_MAX_PATH_LEN] = "initialcondition.vtk";</div>
<div> ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,vtkfilename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);</div><div> ierr = VecView(u,viewer);CHKERRQ(ierr);</div><div> ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr)</div>
<div> }</div><div><br></div><div> .</div><div> .</div><div> ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);</div><div> .</div><div> .</div><div><br></div><div> if (write_output) {</div><div> PetscViewer viewer;</div>
<div> const char *name;</div><div> Vec lv;</div><div> DMGetLocalVector(dm, &lv);</div><div> ierr = PetscObjectGetName((PetscObject) u,&name);CHKERRQ(ierr);</div><div> ierr = PetscObjectSetName((PetscObject) lv, name);CHKERRQ(ierr);</div>
<div> ierr = DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr);</div><div> ierr = DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);CHKERRQ(ierr);</div><div> char vtkfilename[PETSC_MAX_PATH_LEN] = "solution.vtk";</div>
<div> ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,vtkfilename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);</div><div> ierr = VecView(lv,viewer);CHKERRQ(ierr);</div><div> ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); </div>
<div> ierr = DMRestoreLocalVector(dm, &lv);CHKERRQ(ierr);</div><div> }</div><div><br></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Jan 29, 2014 at 2:04 PM, Geoffrey Irving <span dir="ltr"><<a href="mailto:irving@naml.us" target="_blank">irving@naml.us</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>On Wed, Jan 29, 2014 at 12:00 PM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br>
> 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.<br>
><br>
> 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.<br>
<br>
</div>I'll look forward to the improved interface. For better or worse, I'd<br>
like to be able to view simulations in the near future, so it looks<br>
like I have to go with Matt's perverse but working version now.<br>
<span><font color="#888888"><br>
Geoffrey<br>
</font></span><div><div><br>
> On January 29, 2014 12:53:26 PM MST, Geoffrey Irving <<a href="mailto:irving@naml.us" target="_blank">irving@naml.us</a>> wrote:<br>
>>On Wed, Jan 29, 2014 at 11:36 AM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br>
>>> Matt's sample code doesn't set it either. We need to fix (and the<br>
>>*only* acceptable fix if that VecView always does the right thing,<br>
>>because we have to be able to call it in analysis settings that know<br>
>>nothing about your discretization).<br>
>><br>
>>Matt's sample code doesn't set it either, but for Matt's sample code I<br>
>>know where to insert the one line call to DMPlexProjectFunctionLocal.<br>
>>For your version I never have explicit access to the local vector, so<br>
>>I can't insert the fix.<br>
>><br>
>>> The problem is that some vectors reside in a homogeneous space (e.g.<br>
>>increments and eigenvectors) while others reside in the inhomogeneous<br>
>>space (solutions). We can add a flag or BC attribute on the vector to<br>
>>this effect, but this (and slip conditions) was the issue that led me<br>
>>to conclude that removing boundary nodes was mostly a false economy.<br>
>><br>
>>To leave the boundary conditions in, we would need efficient support<br>
>>for a very large, very sparse MatNullSpace. This is doable via<br>
>>shells, but is it easy to do in a way that doesn't interfere with the<br>
>>user's other null spaces?<br>
>><br>
>>Geoffrey<br>
><br>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>