<div dir="ltr"><div><div><div><div>Ok, I found why it was failing. The gmsh reader didn't invert the tetrahedra (or hexahedra) like it does in the exodus reader (both have the same vertex ordering).<br></div>The hdf5 to xdmf issue was because the block size for the local coordinates vector was never set.<br></div>I created a pull request.<br><br></div>Thanks,<br></div>Sander<br></div><div class="gmail_extra"><br><div class="gmail_quote">On 19 June 2016 at 23:03, Sander Arens <span dir="ltr"><<a href="mailto:Sander.Arens@ugent.be" target="_blank">Sander.Arens@ugent.be</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">And the files of course.<br></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On 19 June 2016 at 23:02, Sander Arens <span dir="ltr"><<a href="mailto:Sander.Arens@ugent.be" target="_blank">Sander.Arens@ugent.be</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><div><div><div><div><div><div>Thanks, Matt. That indeed made the problem for the test dissapear.<br><br></div>However, there's still a problem when I try to work with a gmsh mesh.<br></div><br>I modified snes/examples/tutorials/ex77 to read in this mesh (it's just a simple bar generated in gmsh) and I used DMPlexReverseCell because the jacobian determinants where negative. Snes doesn't converge and looking for example in the residual function f1_u_3d shows that x and u are not the same at the beginning (which should be though, as the coordinates where projected onto u). Please see the attached files ex77.c, makefile and bar.msh.<br><br></div>I then tried working with a displacement field instead of deformation and this seemed to work (because the initial values where just zero), but the directions of the face normals are in the wrong direction because the bar was stretched instead of compressed. Please see the attached files ex77_2.c, makefile_2 and bar.msh.<br><br></div>So did I forget something when using a read in mesh or is it something else?<br><br></div>Also, /bin/petsc_gen_xdmf.py seemed to fail on the output ex77_2.h5 of make runex77_2_3 (when using ex77_2.c, makefile_2 and bar.msh).<br><br></div>Thanks for all the help,<br></div>Sander<br></div><div><div><div class="gmail_extra"><br><div class="gmail_quote">On 19 June 2016 at 02:17, 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:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span>On Sat, Jun 18, 2016 at 12:49 PM, Sander Arens <span dir="ltr"><<a href="mailto:Sander.Arens@ugent.be" target="_blank">Sander.Arens@ugent.be</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><span><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote">1) Does the test run for you?</blockquote><br></span>I ran the tests and they all pass.<span><span><br><br></span><span></span><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote">2) If so, then could it be that you assume something about the order of functions?</blockquote></span><div><br>I'm not really sure what you mean with this.<span><br><br><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div>So the values in 'u' are the FEM coefficients for the coordinate function {x, y}. Why are you calling 'x' projected and</div><div>'y' interpolated? </div></blockquote></span></div><span></span></div><div><span><br></span></div><div><span>No, u is a vector with three fields which are obtained in the following way:<br></span></div><div><span><br></span></div><div><span>1) Project the coordinates onto field nr 0 so that I have the coordinates for all the dofs<br></span></div><div><span>2) Project the coordinates onto field nr 1, so that I get the coordinates at the centroid of the cell<br></span></div><div><span>3) Project field nr 0 onto field nr 2, </span><span>so that I get the coordinates at the centroid of the cell<br><br></span></div><div><span>I then output the values of the different fields, so that I can compare the values of field 1 and 2.<br><br></span></div><div><span>What I would expect is that the values of field 1 and 2 would be the same. This is the case for field 0 being P1, but not for P2.<br></span></div><div><span>I also projected the coordinates onto a P2 field in snes/examples/tutorials/ex77.c, but this was done using </span>DMPlexCreateBoxMesh and gave no problems. <br></div><div>Maybe it's some assumption on the vertex numbering in DMPlexCreateFromCellList or DMPlexCreateFromFile that I'm not aware of?<br><br></div><div>I hope I made the problem clearer now?<br></div></div></blockquote><div><br></div></span><div>Okay, I see what happened now. It is indeed a shortcoming. The way I have things setup now, all fields share the quadrature</div><div>since that way I can evaluate all fields in the same loop at a quadrature point. This could of course be changed at some cost.</div><div>I have changed your example to run.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div><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>Thanks,<br></div><div>Sander<br></div><span></span></div><div><div><div class="gmail_extra"><br><div class="gmail_quote">On 18 June 2016 at 16:33, 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:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span>On Sat, Jun 18, 2016 at 9: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-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span>On Sat, Jun 18, 2016 at 8:56 AM, Sander Arens <span dir="ltr"><<a href="mailto:Sander.Arens@ugent.be" target="_blank">Sander.Arens@ugent.be</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div><div><div><div>Hello all,<br><br></div>Recently I've been trying to run a fem problem using DMPlex, PetscFE, etc. using a mesh created with gmsh. I used DMPlexReverseCell to get rid of the negative jacobian determinants. However, I noticed another problem: when I projected the coordinates on a quadratic fem-field and interpolated them they didn't match the interpolated values of the original coordinate field anymore.<br><br></div>I attached a simple test to reproduce this. With linear interpolation I see no problem, it's only when I start using quadratic interpolation.<br><br></div>I have no idea why this is giving problems. Am I missing something simple here or is this a bug?<br></div></div></div></blockquote><div><br></div></span><div>Tests for this already exist</div><div><br></div><div> <a href="https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/src/dm/impls/plex/examples/tests/ex3.c?at=master&fileviewer=file-view-default" target="_blank">https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/src/dm/impls/plex/examples/tests/ex3.c?at=master&fileviewer=file-view-default</a></div><div><br></div><div>with run parameters here:</div><div><br></div><div> <a href="https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/config/builder.py?at=master&fileviewer=file-view-default#builder.py-58" target="_blank">https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/config/builder.py?at=master&fileviewer=file-view-default#builder.py-58</a></div><div><br></div><div>and specifically for P2 triangles</div><div><br></div><div> <a href="https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/config/builder.py?at=master&fileviewer=file-view-default#builder.py-72" target="_blank">https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/config/builder.py?at=master&fileviewer=file-view-default#builder.py-72</a></div><div><br></div><div>1) Does the test run for you?</div><div><br></div><div> ./config/builder2.py check src/dm/impls/plex/examples/tests/ex3.c</div><div><br></div><div>2) If so, then could it be that you assume something about the order of functions?</div></div></div></div></blockquote><div><br></div></span><div>I am having a hard time understanding what you are doing here:</div><div><br></div><div><pre> /* View the projected coordinates field */
ierr = VecGetSubVector(u, fields[0], &subvec);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) subvec, fieldNames[0]);CHKERRQ(ierr);
ierr = VecViewFromOptions(subvec, NULL, "-projectedcoordinates_vec_view");CHKERRQ(ierr);
ierr = VecRestoreSubVector(u, fields[0], &subvec);CHKERRQ(ierr);
/* View the interpolated coordinates field */
ierr = VecGetSubVector(u, fields[1], &subvec);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) subvec, fieldNames[1]);CHKERRQ(ierr);
ierr = VecViewFromOptions(subvec, NULL, "-interpolatedcoordinates_vec_view");CHKERRQ(ierr);
ierr = VecRestoreSubVector(u, fields[1], &subvec);CHKERRQ(ierr);</pre></div><div><br></div><div>So the values in 'u' are the FEM coefficients for the coordinate function {x, y}. Why are you calling 'x' projected and</div><div>'y' interpolated?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><span><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><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-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div></div>Thanks,<br></div>Sander<span><font color="#888888"><br></font></span></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br><br clear="all"><div><br></div>-- <br><div data-smartmail="gmail_signature">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>
</font></span></div></div>
</blockquote></span></div><span><br><br clear="all"><div><br></div>-- <br><div data-smartmail="gmail_signature">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>
</span></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div></div></div><div><div><br><br clear="all"><div><br></div>-- <br><div data-smartmail="gmail_signature">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></div></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>