<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Aug 1, 2014 at 6:01 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">Never mind, I sort of figured it out. What I was kind of looking for was within the DMPlexMatSetClosure(...) function (plex.c). Specifically, these lines and the private functions within:<br>
<br><pre width="80"><a name="14793ce4964e7552_14793b3a9d166246_line5383">5383: </a> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetWorkArray.html#DMGetWorkArray" target="_blank">DMGetWorkArray</a>(dm, numIndices, PETSC_INT, &indices);
<a name="14793ce4964e7552_14793b3a9d166246_line5384">5384: </a> <font color="#4169E1">if</font> (numFields) {
<a name="14793ce4964e7552_14793b3a9d166246_line5385">5385: </a> <font color="#4169E1">for</font> (p = 0; p < numPoints*2; p += 2) {
<a name="14793ce4964e7552_14793b3a9d166246_line5386">5386: </a> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscInt.html#PetscInt" target="_blank">PetscInt</a> o = points[p+1];
<a name="14793ce4964e7552_14793b3a9d166246_line5387">5387: </a> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/IS/PetscSectionGetOffset.html#PetscSectionGetOffset" target="_blank">PetscSectionGetOffset</a>(globalSection, points[p], &globalOff);
<a name="14793ce4964e7552_14793b3a9d166246_line5388">5388: </a> indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PETSC_FALSE.html#PETSC_FALSE" target="_blank">PETSC_FALSE</a>, o, indices);
<a name="14793ce4964e7552_14793b3a9d166246_line5389">5389: </a> }
<a name="14793ce4964e7552_14793b3a9d166246_line5390">5390: </a> } <font color="#4169E1">else</font> {
<a name="14793ce4964e7552_14793b3a9d166246_line5391">5391: </a> <font color="#4169E1">for</font> (p = 0, off = 0; p < numPoints*2; p += 2) {
<a name="14793ce4964e7552_14793b3a9d166246_line5392">5392: </a> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscInt.html#PetscInt" target="_blank">PetscInt</a> o = points[p+1];
<a name="14793ce4964e7552_14793b3a9d166246_line5393">5393: </a> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/IS/PetscSectionGetOffset.html#PetscSectionGetOffset" target="_blank">PetscSectionGetOffset</a>(globalSection, points[p], &globalOff);
<a name="14793ce4964e7552_14793b3a9d166246_line5394">5394: </a> indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PETSC_FALSE.html#PETSC_FALSE" target="_blank">PETSC_FALSE</a>, o, indices);
<a name="14793ce4964e7552_14793b3a9d166246_line5395">5395: </a> }
<a name="14793ce4964e7552_14793b3a9d166246_line5396">5396: </a> }<br><br></pre><pre width="80"><span style="font-family:verdana,sans-serif"><font>If I understand this part of the code correctly, it returns the global indices based on the global offsets; if a certain global dof happens to be constrained the indices will be negative. All I want is the ability to extract global column indices because each row of my equality constraint matrix will be element specific. Is there an existing function that does something like this, or do I have to write it myself?</font></span></pre>
</div></blockquote><div>I think you have to write that one yourself. The closest function is probably</div><div><br></div><div> PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[])</div>
<div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div class="gmail_extra">
<div class="gmail_quote">On Fri, Aug 1, 2014 at 1:40 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">
<div>On Fri, Aug 1, 2014 at 2:38 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div><div><div>As far as I know this equality constraint should exist outside of the large square matrix.<br>
<br></div>What I was thinking of doing is manually creating an <i>Nele</i> by <i>Ndof</i> matrix, each row will have the same number of non-zeros. Since the size of <i>Ndof</i> is simply the global number of free dofs formulated by the PetscSection, I want to know what's the best way to map the free dofs into the global column indices [0 <i>Ndof</i>)<br>
<br></div>For example, in src/dm/impls/plex/examples/tutorials/ex1.c where a box mesh with three fields is created, there are a total of 41 dofs to which 8 of them are constrained. If no constraints were set, then the size of u (i.e. <i>Ndof</i>) would be 41, but with constraints set it is 33. When you view the constrained global vector u, it lists all the dofs in the order as specified by the PetscSection but skips all the constrained ones. Would I have to use something like DMSet/GetDefaultGlobalSection to do the mapping?<br>
</div></div></blockquote><div><br></div></div><div>I do not understand what you are asking. However, you can just use the PetscSection for the Ndof field to give you</div><div>the local sizes for your matrix.</div><div><br>
</div>
<div> Matt</div><div><div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div>
</div>Thanks,<br>Justin<br></div><div class="gmail_extra">
<br><br><div class="gmail_quote">
On Fri, Aug 1, 2014 at 7:26 AM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">
<div>On Wed, Jul 30, 2014 at 5:01 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div><div dir="ltr"><div><div>That's good to know thanks. I have another question somewhat related.<br>
<br>
</div>In that quadratic programming problem I have shown earlier, my K matrix is of global size <i>Ndof</i> by <i>Ndof</i> where <i>Ndof</i>
is number of free dofs (for instance, if I have a 2x2 quadrilateral
mesh and want to evaluate concentration at all the verticies, I will
have nine total dofs but if I have Dirichlet constraints all along the
boundary, only my center vertex is free so <i>Ndof</i> = 1. Thus when I simply call DMCreateMatrix(...) it would give me a 1x1 matrix.)<br>
<br></div>However, I want my M constraint matrix to be of size <i>Nele</i> by <i>Ndof</i> where <i>Nele</i> is the total number of elements. How would I create this matrix? Because ultimately, I am subjecting all of my free u's to <i>Nele</i> number of equality constraints, and if possible I want this M matrix to be compatible with the layout from Vec u.</div>
</div></blockquote><div><br></div></div><div>This is not hard in principle. I used to have code that did this, but I took it out during the rewrite because no one</div><div>was using it and it was complicated. It just generalizes the code in DMPlexPreallocateOperator() to take two</div>
<div>PetscSections instead of one.</div><div><br></div><div>Are you sure that this does not belong inside a large square matrix for the entire problem?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div>
<div>
<div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div class="gmail_extra"><br><div class="gmail_quote">
<div>On Wed, Jul 30, 2014 at 2:03 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
</div><div><div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">
<div><div>On Wed, Jul 30, 2014 at 2:53 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div><div><div><div>Hi,<br><br>This might be a silly question, but I am developing a finite element code to solve the non-negative diffusion equation. I am using the DMPlex structure and plan on subjecting the problem to a bounded constraint optimization to ensure non-negative concentrations. Say for instance, I want to solve:<br>
<br></div>min || K*u - F ||^2<br></div>subject to u >= 0<br></div><div> M*u = 0<br></div><div><br></div>Where K is the Jacobian, F the residual, u the concentration, and M some equality constraint matrix. I know how to do this via Tao, but my question is can Tao handle Mats and Vecs created via DMPlex? Because from the SNES and TS examples I have seen in ex12, ex62, and ex11, it seems there are solvers tailored to handle DMPlex created data structures.<br>
</div></div></blockquote><div><br></div></div></div><div>Yes, TAO can handle objects created by DMPlex since they are just Vec and Mat structures. The additions</div><div>to ex12, ex62, and ex11 concern the callbacks from the solver to the constructions routines for residual and</div>
<div>Jacobian. Right now, none of the callbacks are automatic for TAO so you would have to code them yourself,</div><div>probably by just calling the routines from SNES or TS so it should not be that hard. We are working to get</div>
<div>everything integrated as it is in the other solvers.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div dir="ltr"><div></div>Thanks,<br>Justin<span><font color="#888888"><br></font></span></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br><br clear="all"><span><font color="#888888"><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
</font></span></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><br></font></span></div><span><font color="#888888">
</font></span></blockquote></div></div><span><font color="#888888"><div><br><br clear="all"><span class=""><font color="#888888"><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
</font></span></div></font></span></div></div><span class=""><font color="#888888">
</font></span></blockquote></div><span class=""><font color="#888888"><br></font></span></div><span class=""><font color="#888888">
</font></span></blockquote></div></div></div><span class=""><font color="#888888"><div><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></font></span></div></div>
</blockquote></div><br></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>