# [petsc-users] DMPLEX project function error

Matthew Knepley knepley at gmail.com
Tue Dec 18 08:33:33 CST 2018

On Sun, Oct 14, 2018 at 12:36 AM Ellen M. Price <ellen.price at cfa.harvard.edu>
wrote:

> Hi Matt,
>
> I ran into a problem implementing this. My integration code relies on
> knowing the points of each tetrahedron in the domain so I can turn an
> integral over all space into a sum of integrals over cells. I can't find
> an example of getting vertices for a given cell (I think TS ex11 comes
> closest, looks like it gets faces?). Is there an elegant way to do this,
> or an example I can follow?
>

I do this in the library, but you can do it yourself. I assume what you
mean is getting
the coordinates of the vertices of the tet. If you only have coordinates at
vertices, then
its easy. You just use

https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexVecGetClosure.html

to get the coordinates over the cell. We use this in all the geometry
routines.

Thanks,

Matt

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

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their