[petsc-users] DMPLEX project function error
Ellen M. Price
ellen.price at cfa.harvard.edu
Sat Oct 13 23:36:19 CDT 2018
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?
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
> their experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/%7Eknepley/>
More information about the petsc-users
mailing list