<div dir="ltr"><div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Sun, Oct 14, 2018 at 12:36 AM Ellen M. Price <<a href="mailto:ellen.price@cfa.harvard.edu">ellen.price@cfa.harvard.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi Matt,<br>
<br>
I ran into a problem implementing this. My integration code relies on<br>
knowing the points of each tetrahedron in the domain so I can turn an<br>
integral over all space into a sum of integrals over cells. I can't find<br>
an example of getting vertices for a given cell (I think TS ex11 comes<br>
closest, looks like it gets faces?). Is there an elegant way to do this,<br>
or an example I can follow?<br></blockquote><div><br></div><div>I do this in the library, but you can do it yourself. I assume what you mean is getting</div><div>the coordinates of the vertices of the tet. If you only have coordinates at vertices, then</div><div>its easy. You just use</div><div><br></div><div> <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexVecGetClosure.html">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexVecGetClosure.html</a></div><div><br></div><div>to get the coordinates over the cell. We use this in all the geometry routines.</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">
Ellen<br>
<br>
<br>
On 10/13/2018 04:26 AM, Matthew Knepley wrote:<br>
> On Fri, Oct 12, 2018 at 6:00 PM Ellen M. Price<br>
> <<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>> wrote:<br>
> <br>
> Hi Matt,<br>
> <br>
> Thanks for the feedback! I'm always open to editorial comments. Do you<br>
> have a suggestion for a better way? I'm basically just trying to do an<br>
> integral over all space to get my boundary values for the FEM solve. I<br>
> don't know if PETSc has a way to do something like this?<br>
> <br>
> <br>
> So you are defining your BC by a boundary integral. I would just stick<br>
> that integration code<br>
> into a function and let that be called as your BC function. It will not<br>
> be called "extra" times,<br>
> just once for each dual vector, which in the case of P1 is every vertex.<br>
> This way it will work<br>
> for any element.<br>
> <br>
> Thanks,<br>
> <br>
> Matt<br>
> <br>
> <br>
> Ellen<br>
> <br>
> <br>
> On 10/12/2018 04:36 PM, Matthew Knepley wrote:<br>
> > On Fri, Oct 12, 2018 at 3:51 PM Ellen M. Price<br>
> > <<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>>> wrote:<br>
> ><br>
> > Hi Matt,<br>
> ><br>
> > Basically, I have a simple meshing script that generates a set<br>
> of points<br>
> > for my FEM mesh and outputs a C file that compiles into my<br>
> program.<br>
> > Along with those points is a static data field that I need to be<br>
> > associated with the correct data point. I don't have a<br>
> function that<br>
> > evaluates at each point, it's just pre-generated data.<br>
> ><br>
> ><br>
> > I cannot help but make editorial comments. Do not do this. It is<br>
> > incredibly fragile, since<br>
> > it depends on both a mesh and a discretization.<br>
> > <br>
> ><br>
> > Since I called DMPlexDistribute() early on, I think my mesh<br>
> points get<br>
> > re-ordered, right? So I need a way to re-order my data field<br>
> to match<br>
> > that ordering. Is there a good way to do this?<br>
> ><br>
> ><br>
> > However, if you want to do this, it is not hard. You do use<br>
> > DMPlexDistributeField<br>
> ><br>
> (<a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexDistributeField.html" rel="noreferrer" target="_blank">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexDistributeField.html</a>),<br>
> > but you need<br>
> > to be careful about the arguments. You want<br>
> ><br>
> > DMPlexDistribute(dmOriginal, 0, &migrationSF, &dmNew);<br>
> > DMPlexDistributeField(dmOriginal, migrationSF, originalSection,<br>
> > originalVec, newSection, newVec);<br>
> ><br>
> > Your Sections are easy to make, just put 1 dof on every vertex. The<br>
> > originalVec is just a Vec of your values<br>
> > (you can use VecCreateSeqWithArray() to wrap one around it), and the<br>
> > newVec you can get with<br>
> > DMCreateLocalVector(dmNew, &newVec) if you do DMSetSection(dmNew,<br>
> > newSection).<br>
> ><br>
> > Does this make sense?<br>
> ><br>
> > Thanks,<br>
> ><br>
> > Matt<br>
> > <br>
> ><br>
> > I tried using DMPlexDistributeField as a first attempt, but I<br>
> get some<br>
> > kind of error no matter what I do (presumably I'm not giving<br>
> the right<br>
> > input), and there aren't any examples in the docs.<br>
> ><br>
> > Ellen<br>
> ><br>
> ><br>
> > On 10/12/2018 11:31 AM, Matthew Knepley wrote:<br>
> > > On Fri, Oct 12, 2018 at 11:12 AM Ellen M. Price<br>
> > > <<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>><br>
> > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>><br>
> > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>>>> wrote:<br>
> > ><br>
> > > So I found *that* problem -- turns out that using<br>
> > DMPlexCreateSection<br>
> > > was messing things up, and it seems PETSc is handling<br>
> that part<br>
> > > internally?<br>
> > ><br>
> > ><br>
> > > DMPlexSection creates the Section based on the discretization<br>
> > > information in the DM (stored in a PetscDS).<br>
> > > If no discretization info is available, then it defaults to P0 I<br>
> > think,<br>
> > > which is what you got.<br>
> > ><br>
> > > After discretization info is put into the DM, if I ask for the<br>
> > Section,<br>
> > > the DM will create it automatically if it is missing.<br>
> > > <br>
> > ><br>
> > > However, I've run into a related problem now. I need to<br>
> > > project some field data onto the mesh; it's a static<br>
> array of<br>
> > > integration weights, and each one corresponds to a single<br>
> > point in the<br>
> > > mesh. Based on this answer:<br>
> > ><br>
> > > <br>
> > <br>
> <a href="https://lists.mcs.anl.gov/pipermail/petsc-users/2015-December/027964.html" rel="noreferrer" target="_blank">https://lists.mcs.anl.gov/pipermail/petsc-users/2015-December/027964.html</a><br>
> > ><br>
> > > I think I need DMPlexDistributeField.<br>
> > ><br>
> > ><br>
> > > I don't think that is what you want. That is just sending data<br>
> > around.<br>
> > ><br>
> > > You want to define an FEM field. We use DMProject*() for<br>
> this. These<br>
> > > take as input either<br>
> > > a function of the coordinates (regular function) or another<br>
> FEM field.<br>
> > ><br>
> > > Can you tell me how your function is defined? I mean apart<br>
> from the<br>
> > > particular data you<br>
> > > have at points.<br>
> > ><br>
> > > Thanks,<br>
> > ><br>
> > > Matt<br>
> > > <br>
> > ><br>
> > > However, *that* function takes a<br>
> > > section as an argument, and I had to take out the part<br>
> where I<br>
> > create a<br>
> > > section to get things working! Calling<br>
> DMGetDefaultSection or<br>
> > > DMGetDefaultGlobalSection doesn't seem to help matters,<br>
> > either, because<br>
> > > then I seem to get a null section?<br>
> > ><br>
> > > So I either need to figure out how to set the section<br>
> properly<br>
> > or how to<br>
> > > get the section properly. As always, any help is<br>
> appreciated.<br>
> > ><br>
> > > Ellen<br>
> > ><br>
> > ><br>
> > > On 10/11/2018 05:41 PM, Matthew Knepley wrote:<br>
> > > > On Thu, Oct 11, 2018 at 4:39 PM Ellen M. Price<br>
> > > > <<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>><br>
> > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>><br>
> > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>><br>
> > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>>><br>
> > > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>><br>
> > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>><br>
> > > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>><br>
> > <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a><br>
> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>>>>> wrote:<br>
> > > ><br>
> > > > I was working with a DMPLEX and FEM, following<br>
> SNES example<br>
> > > 12. I get<br>
> > > > the following error when I call DMProjectFunction,<br>
> but I<br>
> > don't<br>
> > > know what<br>
> > > > it means. Can anyone explain where I might have gone<br>
> > wrong, or<br>
> > > at least<br>
> > > > what this error is telling me? I think the point<br>
> closure<br>
> > size is<br>
> > > > correct, since my mesh is 3d simplex,<br>
> > > ><br>
> > > ><br>
> > > > Yes, if you have 3D SIMPLEX mesh and are using P1<br>
> elements,<br>
> > then you<br>
> > > > would have<br>
> > > > 4 dofs in the closure of a cell. The dual space<br>
> dimension is the<br>
> > > number<br>
> > > > of dual space<br>
> > > > basis vectors assigned to points in the closure. Since<br>
> it is<br>
> > 1, it<br>
> > > looks<br>
> > > > like you have a P0<br>
> > > > dual space. I assume you changed something in ex12?<br>
> > > ><br>
> > > > Thanks,<br>
> > > ><br>
> > > > Matt<br>
> > > > <br>
> > > ><br>
> > > > but what is the dual space<br>
> > > > dimension, and where might I have set it incorrectly?<br>
> > > ><br>
> > > > [0]PETSC ERROR: Nonconforming object sizes<br>
> > > > [0]PETSC ERROR: The section point closure size 4<br>
> != dual<br>
> > space<br>
> > > > dimension 1<br>
> > > > [0]PETSC ERROR: See<br>
> > > <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a><br>
> > > > for trouble shooting.<br>
> > > > [0]PETSC ERROR: Petsc Release Version 3.9.2, May,<br>
> 20, 2018<br>
> > > > ...<br>
> > > > [0]PETSC ERROR: Configure options<br>
> > > > --prefix=/home/eprice/software/petsc-opt --with-hdf5=1<br>
> > > > --with-hdf5-dir=/home/eprice/software/hdf5-parallel<br>
> > --with-mpe=1<br>
> > > > --with-mpe-dir=/home/eprice/software/mpe<br>
> --with-debugging=0<br>
> > > > LDFLAGS="-pthread -lz" COPTFLAGS="-O3 -march=native<br>
> > -mtune=native"<br>
> > > > CXXOPTFLAGS="-O3 -march=native -mtune=native"<br>
> FOPTFLAGS="-O3<br>
> > > > -march=native -mtune=native" --with-mpi=1<br>
> > > > --with-mpi-dir=/home/eprice/software/mpich<br>
> --with-mumps=1<br>
> > > > --with-mumps-dir=/home/eprice/software/mumps<br>
> > --with-parmetis=1<br>
> > > > --with-parmetis-dir=/home/eprice/software/parmetis<br>
> > --with-metis=1<br>
> > > > --with-metis-dir=/home/eprice/software/parmetis<br>
> > --with-ptscotch=1<br>
> > > > --with-ptscotch-dir=/home/eprice/software/scotch<br>
> > > --with-scalapack=1<br>
> > > > --with-scalapack-dir=/home/eprice/software/scalapack<br>
> > > > [0]PETSC ERROR: #1 DMProjectLocal_Generic_Plex()<br>
> line 347 in<br>
> > > > <br>
> /h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c<br>
> > > > [0]PETSC ERROR: #2 DMProjectFunctionLocal_Plex()<br>
> line 428 in<br>
> > > > <br>
> /h/sabriel0/src/petsc-3.9.2/src/dm/impls/plex/plexproject.c<br>
> > > > [0]PETSC ERROR: #3 DMProjectFunctionLocal() line<br>
> 6265 in<br>
> > > > /h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c<br>
> > > > [0]PETSC ERROR: #4 DMProjectFunction() line 6250 in<br>
> > > > /h/sabriel0/src/petsc-3.9.2/src/dm/interface/dm.c<br>
> > > > ...<br>
> > > ><br>
> > > > (I know this is an optimized PETSc build, but I<br>
> get the same<br>
> > > error from<br>
> > > > my debug build, it's just much slower.)<br>
> > > ><br>
> > > > Ellen<br>
> > > ><br>
> > > ><br>
> > > ><br>
> > > > --<br>
> > > > What most experimenters take for granted before they<br>
> begin their<br>
> > > > experiments is infinitely more interesting than any<br>
> results<br>
> > to which<br>
> > > > their experiments lead.<br>
> > > > -- Norbert Wiener<br>
> > > ><br>
> > > > <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> <<a href="https://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> > <<a href="https://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> > > <<a href="https://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> > > > <<a href="http://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> > ><br>
> > ><br>
> > ><br>
> > > --<br>
> > > What most experimenters take for granted before they begin their<br>
> > > experiments is infinitely more interesting than any results<br>
> to which<br>
> > > their experiments lead.<br>
> > > -- Norbert Wiener<br>
> > ><br>
> > > <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> <<a href="https://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> > <<a href="https://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> > > <<a href="http://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> ><br>
> ><br>
> ><br>
> > --<br>
> > What most experimenters take for granted before they begin their<br>
> > experiments is infinitely more interesting than any results to which<br>
> > their experiments lead.<br>
> > -- Norbert Wiener<br>
> ><br>
> > <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> <<a href="https://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> > <<a href="http://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/%7Eknepley/</a>><br>
> <br>
> <br>
> <br>
> -- <br>
> What most experimenters take for granted before they begin their<br>
> experiments is infinitely more interesting than any results to which<br>
> their experiments lead.<br>
> -- Norbert Wiener<br>
> <br>
> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> <<a href="http://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/%7Eknepley/</a>><br>
</blockquote></div><br clear="all"><div><br></div>-- <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="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div></div>