<div dir="ltr"><div dir="ltr">On Fri, Jun 20, 2025 at 5:14 PM Noam T. <<a href="mailto:dontbugthedevs@proton.me">dontbugthedevs@proton.me</a>> wrote:</div><div class="gmail_quote gmail_quote_container"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px">Thank you once again, the code provides exactly what needed. <br></div><div style="font-family:Arial,sans-serif;font-size:14px">An alternative for the VTK use was subdividing cells using corner nodes and integration points, such that all cells were first order. Any "better" alternative format/visualization software for this purpose?</div></blockquote><div><br></div><div>I do have this implemented. If you give -dm_plex_high_order_view, it will refine the grid and project into the linear space.</div><div>You can see that the code is pretty simple:</div><div><br></div><div>  <a href="https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/plex.c?ref_type=heads*L2021__;Iw!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-wmRROXs$">https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/plex.c?ref_type=heads#L2021</a></div><div><br></div><div>It uses DMPlexInterpolate() to connect the spaces.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px">Possibly the last question in relation to this matter. We use first order meshes only, as you suggest, and let PETSc handle everything high-order through the approximation space. Hence, when retrieving node coordinates with DMGetCoordinates(Local) or DMPlexGetCellCoordinates, one gets the corner nodes only. </div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">Is the list of additional, high-order nodes coordinates readily available (stored) somewhere to be retrieved?<br></div></blockquote><div><br></div><div>Yes. The idea is to think of coordinates as a discretized field on the mesh, exactly as the solution field. Thus if you want higher order coordinates, you choose a higher order coordinate space. I give the coordinate space prefix cdm_, so you could say</div><div><br></div><div>  -cdm_petscspace_degree 2</div><div><br></div><div>to get quadratic coordinates. There are some tests in Plex tests ex33.c</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:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px"></div><div style="font-family:Arial,sans-serif;font-size:14px">They can be computed (e.g. using DMPlexReferenceToCoordinates knowing their position in the reference cell; or using corner nodes coordinates), but this will result in shared nodes being computed possibly several times; the large the mesh, the worse.</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div>E.g. in the example of the image attached  before, DMPlexGetClosureIndices returns </div><div>[4, 5, 6, 0, 1, 3] for the first cell</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">[7, 8, 5, 1, 2, 3] for the second cell<br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">so that 3 nodes (5, 1, 3) will be computed twice if done naively, cell by cell.</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Thank you,</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Noam<br></div><div>
        On Thursday, June 19th, 2025 at 12:43 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
        <blockquote type="cite">
            <div dir="ltr"><div dir="ltr">On Wed, Jun 18, 2025 at 6:49 PM Noam T. <<a href="mailto:dontbugthedevs@proton.me" rel="noreferrer nofollow noopener" target="_blank">dontbugthedevs@proton.me</a>> wrote:</div><div class="gmail_quote"><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div>See image attached.</div><div>Connectivity of the top mesh (first order triangle), can be obtained with the code shared before.</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Connectivity of the  bottom mesh (second order triangle) is what I would be interested in obtaining.<br></div><div><span style="font-family:Arial,sans-serif"><br></span></div><div><span style="font-family:Arial,sans-serif">However, given your clarification on what the Plex and the PetscSection handle, it might not work; I am trying to get form the Plex what's only available from the PetscSection.<br></span></div><div><span style="font-family:Arial,sans-serif"><br></span></div><div><span style="font-family:Arial,sans-serif">The purpose of this extended connectivity is plotting; in particular, using VTU files, where the "connectivity" of cells is required, and the extra nodes would be needed when using higher-order elements (e.g. VTK_QUADRATIC_TRIANGLE, VTK_QUADRATIC_QUAD, etc).</span></div></blockquote><div><br></div><div>Oh yes. VTK does this in a particularly ugly and backward way. Sigh. There is nothing we can do about this now, but someone should replace VTK with a proper interface at some point.</div><div><br></div><div>So I understand why you want it and it is a defensible case, so here is how you get that (with some explanation). Those locations, I think, should not be understood as topological things, but rather as the locations of point evaluation functionals constituting a basis for the dual space (to your approximation space). I would call DMPlexGetClosureIndices() (<a href="https://urldefense.us/v3/__https://petsc.org/main/manualpages/DMPlex/DMPlexGetClosureIndices/__;!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0p5NBkz$" rel="noreferrer nofollow noopener" target="_blank">https://petsc.org/main/manualpages/DMPlex/DMPlexGetClosureIndices/</a>) with a Section having the layout of P2 or Q2. This is the easy way to make that</div><div><br></div><div>  PetscSection gs;</div>  PetscFE        fe;<br>  DMPolytopeType ct;<br>  PetscInt       dim, cStart;<br><div><br></div><div>  PetscCall(DMGetDimension(dm, &dim));<br>  PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));<br>  PetscCall(DMPlexGetCellType(dm, cStart, &ct));<br>  PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 2, PETSC_DETERMINE, &fe));<br>  PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));<br>  PetscCall(PetscFEDestroy(&fe));<br></div><div>  PetscCall(DMCreateDS(dm));<br></div><div>  PetscCall(DMGetGlobalSection(dm, &gs));</div><div><br></div><div>  PetscInt *indices = NULL;</div><div>  PetscInt  Nidx;</div><div><br></div><div>  PetscCall(DMPlexGetClosureIndices(dm, gs, gs, cell, PETSC_TRUE, &Nidx, &indices, NULL, NULL));</div><div><br></div><div>  Thanks,</div><div><br></div><div>     MAtt</div><div> </div><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div><span style="font-family:Arial,sans-serif">Perhaps I am over complicating things, and all this information can be obtained in a different, simpler way.<br></span></div><div><span></span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace"><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace">Thanks.</span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace">Noam<br></span></div><div>
        On Tuesday, June 17th, 2025 at 5:42 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" rel="noreferrer nofollow noopener" target="_blank">knepley@gmail.com</a>> wrote:<br>
        <blockquote type="cite">
            <div dir="ltr"><div dir="ltr"><div dir="ltr">On Tue, Jun 17, 2025 at 12:43 PM Noam T. <<a rel="noreferrer nofollow noopener" href="mailto:dontbugthedevs@proton.me" target="_blank">dontbugthedevs@proton.me</a>> wrote:</div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px">Thank you.  For now, I am dealing with vertices only.</div><div style="font-family:Arial,sans-serif;font-size:14px"><br>Perhaps I did not explain myself properly, or I misunderstood your response.</div><div style="font-family:Arial,sans-serif;font-size:14px">What I meant to say is, given an element of order higher than one, the connectivity matrix I obtain this way only contains as many entries as the first order element: 3 for a triangle, 4 for a tetrahedron, etc.</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">Looking at the closure of any cell in the mesh, this is also the case.However, the nodes are definitely present; e.g. from</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px"><span>DMPlexGetCellCoordinates</span>(dm, cell, NULL, nc, NULL, NULL)</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">nc returns the expected value (12 for a 2nd order 6-node planar triangle, 30 for a 2nd order 10-node tetrahedron, etc).</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">The question is, are the indices of these extra nodes obtainable in a similar way as with the code shared before? So that one can have e.g.  [0, 1, 2, 3, 4, 5] for a second order triangle, not just [0, 1, 2].</div></blockquote><div><br></div><div>I am having a hard time understanding what you are after. I think this is because many FEM approaches confuse topology with analysis.</div><div><br></div><div>The Plex stores topology, and you can retrieve adjacencies between any two mesh points.</div><div><br></div><div>The PetscSection maps mesh points (cells, faces, edges , vertices) to sets of dofs. This is how higher order elements are implemented. Thus, we do not have to change topology to get different function spaces.</div><div><br></div><div>The intended interface is for you to call DMPlexVecGetClosure() to get the closure of a cell (or face, or edge). You can also call DMPlexGetClosureIndices(), but index wrangling is what I intended to eliminate.</div><div><br></div><div>What exactly are you looking for here?</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:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px">Thank you.</div><div style="font-family:Arial,sans-serif;font-size:14px">Noam<br></div><div>
        On Friday, June 13th, 2025 at 3:05 PM, Matthew Knepley <<a rel="noreferrer nofollow noopener" href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
        <blockquote type="cite">
            <div dir="ltr"><div dir="ltr">On Thu, Jun 12, 2025 at 4:26 PM Noam T. <<a href="mailto:dontbugthedevs@proton.me" rel="noreferrer nofollow noopener" target="_blank">dontbugthedevs@proton.me</a>> wrote:</div><div class="gmail_quote"><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div><br></div><div>Thank you for the code; it provides exactly what I was looking for.</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Following up on this matter, does this method not work for higher order elements? For example, using an 8-node quadrilateral, exporting to a <span>PETSC_VIEWER_HDF5_VIZ</span> viewer provides the correct matrix of node coordinates in geometry/vertices</div></blockquote><div><br></div><div>If you wanted to include edges/faces, you could do it. First, you would need to decide how you would number things For example, would you number all points contiguously, or separately number cells, vertices, faces and edges. Second, you would check for faces/edges in the closure loop. Right now, we only check for vertices.</div><div><br></div><div>I would say that this is what convinced me not to do FEM this way.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">(here a quadrilateral in [0, 10])</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span>5.0,   5.0</span><div><span>0.0,     0.0</span></div><div><span>10.0,        0.0</span></div><div><span>10.0,        10.0</span></div><div><span>0.0,        10.0</span></div><div><span>5.0, 0.0</span></div><div><span>10.0,       5.0</span></div><div><span>5.0, 10.0</span></div><div><span>0.0,        5.0</span></div><span></span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">but the connectivity in viz/topology is</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">0 1 2 3</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">which are likely the corner nodes of the initial, first-order element, before adding extra nodes for the higher degree element.</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">This connectivity values [0, 1, 2, 3, ...] are always the same, including for other elements, whereas the coordinates are correct</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">E.g. for 3rd order triangle in [0, 1], coordinates are given left to right, bottom to top<br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">0, 0</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">1/3, 0,</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">2/3, 0,</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">1, 0</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">0, 1/3</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">1/3, 1/3</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">2/3, 1/3</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">0, 2/3,</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">1/3, 2/3</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">0, 1<br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">but the connectivity (viz/topology/cells) is  [0, 1, 2].</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Test meshes were created with gmsh from the python API, using</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span>gmsh.option.setNumber("Mesh.ElementOrder", n), for n = 1, 2, 3, ...</span><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Thank you.</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Noam<br></div><div>
        On Friday, May 23rd, 2025 at 12:56 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com" rel="noreferrer nofollow noopener" target="_blank">knepley@gmail.com</a>> wrote:<br>
        <blockquote type="cite">
            <div dir="ltr"><div dir="ltr">On Thu, May 22, 2025 at 12:25 PM Noam T. <<a rel="noreferrer nofollow noopener" href="mailto:dontbugthedevs@proton.me" target="_blank">dontbugthedevs@proton.me</a>> wrote:</div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px">Hello,</div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">Thank you the various options. <br></div><div style="font-family:Arial,sans-serif;font-size:14px"><br></div><div style="font-family:Arial,sans-serif;font-size:14px">Use case here would be obtaining the exact output generated by option 1), DMView() with PETSC_VIEWER_HDF5_VIZ; in particular, the matrix generated under /viz/topology/cells. </div><div style="font-family:Arial,sans-serif;font-size:14px"><br>
        <blockquote type="cite">
            <div dir="ltr"><div><div>There are several ways you might do this. It helps to know what you are aiming for.</div><div><br></div><div>1)
 If you just want this output, it might be easier to just DMView() with
the PETSC_VIEWER_HDF5_VIZ format, since that just outputs the
cell-vertex topology and coordinates</div></div></div></blockquote><div><br></div><div>Is it possible to get this information in memory, onto a Mat, Vec or
some other Int array object directly? it would be handy to have it in order to manipulate it and/or save it to a different format/file. Saving to an HDF5 and loading it
again seems redundant. <br></div><div><br></div><blockquote type="cite"><div dir="ltr"><div><br></div><div><div>2) You can call DMPlexUninterpolate() to produce a mesh with just cells and vertices, and output it in any format.</div><div><br></div><div>3)
 If you want it in memory, but still with global indices (I don't
understand this use case), then you can use DMPlexCreatePointNumbering()
 for an overall global numbering, or DMPlexCreateCellNumbering() and
DMPlexCreateVertexNumbering() for separate global numberings.</div><div><br></div></div></div>

        </blockquote><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Perhaps I missed it, but getting the connectivity matrix in /viz/topology/cells/ did not seem directly trivial to me from the list of global indices returned by  DMPlexGetCell/Point/VertexNumbering() (i.e. I assume all the operations done when calling DMView()).</div></blockquote><div><br></div><div>Something like</div><div><br></div><div>DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);</div><div>DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);</div><div>DMPlexGetVertexNumbering(dm, &globalVertexNumbers);</div><div>ISGetIndices(globalVertexNumbers, &gv);</div><div>for (PetscInt c = cStart; c < cEnd; ++c) {</div><div>    PetscInt *closure = NULL;</div><div><br></div><div>    DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);</div><div>    for (PetscInt cl = 0; c < Ncl * 2; cl += 2) {</div><div>        if (closure[cl] < vStart || closure[cl] >= vEnd) continue;</div><div>        const PetscInt v = gv[closure[cl]] < 0 ? -(gv[closure[cl]] + 1) : gv[closure[cl]];</div><div><br></div><div>        // Do something with v</div><div>    }</div><div>    DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);</div><div>}</div><div>ISRestoreIndices(globalVertexNumbers, &gv);</div><div>ISDestroy(&globalVertexNumbers);</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Thanks,<br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Noam.<br></div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a rel="noreferrer nofollow noopener" href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>

        </blockquote><br>
    </div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div class="gmail_signature" dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$" rel="noreferrer nofollow noopener" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>

        </blockquote><br>
    </div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a rel="noreferrer nofollow noopener" href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</div>

        </blockquote><br>
    </div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div class="gmail_signature" dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$" rel="noreferrer nofollow noopener" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>

        </blockquote><br>
    </div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>