[petsc-users] DMPlex tetrahedra facets orientation

Matthew Knepley knepley at gmail.com
Mon Mar 8 12:22:39 CST 2021


On Mon, Mar 8, 2021 at 11:18 AM Nicolas Barral <
nicolas.barral at math.u-bordeaux.fr> wrote:

> On 08/03/2021 15:55, Matthew Knepley wrote:
> > On Mon, Mar 8, 2021 at 4:02 AM Nicolas Barral
> > <nicolas.barral at math.u-bordeaux.fr
> > <mailto:nicolas.barral at math.u-bordeaux.fr>> wrote:
> >
> >     On 07/03/2021 22:56, Matthew Knepley wrote:
> >      > On Sun, Mar 7, 2021 at 4:51 PM Nicolas Barral
> >      > <nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      > <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>> wrote:
> >      >
> >      >
> >      >     On 07/03/2021 22:30, Matthew Knepley wrote:
> >      >      > On Sun, Mar 7, 2021 at 4:13 PM Nicolas Barral
> >      >      > <nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>
> >      >      > <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>>> wrote:
> >      >      >
> >      >      >     On 07/03/2021 16:54, Matthew Knepley wrote:
> >      >      >      > On Sun, Mar 7, 2021 at 8:52 AM Nicolas Barral
> >      >      >      > <nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>
> >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>>
> >      >      >      > <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>
> >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>>>> wrote:
> >      >      >      >
> >      >      >      >     Matt,
> >      >      >      >
> >      >      >      >     Thanks for your answer.
> >      >      >      >
> >      >      >      >     However, DMPlexComputeCellGeometryFVM does not
> >     compute
> >      >     what I
> >      >      >     need
> >      >      >      >     (normals of height 1 entities). I can't find any
> >      >     function doing
> >      >      >      >     that, is
> >      >      >      >     there one ?
> >      >      >      >
> >      >      >      >
> >      >      >      > The normal[] in DMPlexComputeCellGeometryFVM() is
> >     exactly what
> >      >      >     you want.
> >      >      >      > What does not look right to you?
> >      >      >
> >      >      >
> >      >      >     So it turns out it's not what I want because I need
> >      >     non-normalized
> >      >      >     normals. It doesn't seem like I can easily retrieve
> >     the norm,
> >      >     can I?
> >      >      >
> >      >      >
> >      >      > You just want area-weighted normals I think, which means
> >     that you
> >      >     just
> >      >      > multiply by the area,
> >      >      > which comes back in the same function.
> >      >      >
> >      >
> >      >     Ah by the area times 2, of course, my bad.
> >      >     Do you order height-1 elements in a certain way ? I need to
> >     access the
> >      >     facet (resp. edge) opposite to a vertex in a tet (resp.
> >     triangle).
> >      >
> >      >
> >      > Yes. Now that I have pretty much settled on it, I will put it in
> the
> >      > manual. It is currently here:
> >      >
> >      >
> >
> https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/plexinterpolate.c#L56
> >      >
> >      > All normals are outward facing, but hopefully the ordering in the
> >     sourse
> >      > file makes sense.
> >
> >     Thanks Matt, but I'm not sure I understand well. What I do so far is:
> >
> >     ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
> >         for (i=0; i<dim+1; ++i) {
> >           f = cone[i];
> >           ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, NULL,
> >     &vn[i*dim]);CHKERRQ(ierr);
> >           if (dim == 3) { area *= 2; }
> >           for (j=0; j<dim; ++j) { vn[i*dim+j] *= area; }
> >
> >     So in 3D, it seems that:
> >     (vn[9],vn[10],vn[11]) is the inward normal to the facet opposing
> vertex0
> >     (vn[6],vn[7],vn[8])             "                    "
>    1
> >     (vn[3],vn[4],vn[5])             "                    "
>    2
> >     (vn[0],vn[1],vn[2])             "                    "
>    3
> >
> >     in 2D:
> >     (vn[2],vn[3]) is a normal to the edge opposing vertex 0
> >     (vn[4],vn[5])          "                  "           1
> >     (vn[0],vn[1])          "                  "           2
> >     Yet in 2D, whether the normals are inward or outward does not seem
> >     consistent across elements.
> >
> >     What am I wrongly assuming ?
> >
> >
> > Ah, I see the problem. I probably need another function. You can tell
> > that not many people use Plex this way yet.
> > The logic for what you want is embedded my traversal, but it simple:
> >
> > ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr);
> > ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
> > ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr);
> >     for (i=0; i<coneSize; ++i) {
> >       f = cone[i];
> >       flip = ornt[i] >= 0 ? 1 : -1;
> >       ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, NULL,
> > &vn[i*dim]);CHKERRQ(ierr);
> >       if (dim == 3) { area *= 2; }
> >       for (j=0; j<dim; ++j) { vn[i*dim+j] *= area * flip; }
> > I could make a function that returns all normals, properly oriented. It
> > would just do this.
>
> Ah this works now, thanks Matt. Toby is correct, it is ultimately
> related to Jacobians, and what I need can be done differently, not sure
> it's clearer though.
>
> Out of curiosity, what is the logic in the facet ordering ?
>

The order of faces in a cell was somewhat arbitrary. However, I wanted that

  select vertices from closure(cell) = vertices before interpolation of cell

so the canonical orientation of face should have the vertices such that
they give
me the order of vertices I expect in cell-vertex meshes. This way

  Uninterpolate(Interpolate(dm))

is idempotent.

  Thanks,

     Matt


> Thanks
>
> --
> Nicolas
>
> >
> >    Thanks,
> >
> >       Matt
> >
> >     Thanks,
> >
> >     --
> >     Nicolas
> >
> >      >
> >      >    Thanks,
> >      >
> >      >      Matt
> >      >
> >      >     Thanks
> >      >
> >      >     --
> >      >     Nicolas
> >      >
> >      >
> >      >      >    Thanks,
> >      >      >
> >      >      >      Matt
> >      >      >
> >      >      >     If not, I'll fallback to computing them by hand for
> >     now. Is the
> >      >      >     following assumption safe or do I have to use
> >      >     DMPlexGetOrientedFace?
> >      >      >       >  if I call P0P1P2P3 a tet and note x the cross
> >     product,
> >      >      >       >  P3P2xP3P1 is the outward normal to face P1P2P3
> >      >      >       >  P0P2xP0P3              "                P0P2P3
> >      >      >       >  P3P1xP3P0              "                P0P1P3
> >      >      >       >  P0P1xP0P2              "                P0P1P2
> >      >      >
> >      >      >     Thanks
> >      >      >
> >      >      >     --
> >      >      >     Nicolas
> >      >      >      >
> >      >      >      >    Thanks,
> >      >      >      >
> >      >      >      >      Matt
> >      >      >      >
> >      >      >      >     So far I've been doing it by hand, and after a
> >     lot of
> >      >      >     experimenting the
> >      >      >      >     past weeks, it seems that if I call P0P1P2P3 a
> >     tetrahedron
> >      >      >     and note x
> >      >      >      >     the cross product,
> >      >      >      >     P3P2xP3P1 is the outward normal to face P1P2P3
> >      >      >      >     P0P2xP0P3              "                P0P2P3
> >      >      >      >     P3P1xP3P0              "                P0P1P3
> >      >      >      >     P0P1xP0P2              "                P0P1P2
> >      >      >      >     Have I been lucky but can't expect it to be
> true ?
> >      >      >      >
> >      >      >      >     (Alternatively, there is a link between the
> normals
> >      >     and the
> >      >      >     element
> >      >      >      >     Jacobian, but I don't know the formula and can
> >     find them)
> >      >      >      >
> >      >      >      >
> >      >      >      >     Thanks,
> >      >      >      >
> >      >      >      >     --
> >      >      >      >     Nicolas
> >      >      >      >
> >      >      >      >     On 08/02/2021 15:19, Matthew Knepley wrote:
> >      >      >      >      > On Mon, Feb 8, 2021 at 6:01 AM Nicolas Barral
> >      >      >      >      > <nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>
> >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>>
> >      >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>
> >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>>>
> >      >      >      >      > <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>
> >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>>
> >      >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>
> >      >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>
> >      >     <mailto:nicolas.barral at math.u-bordeaux.fr
> >     <mailto:nicolas.barral at math.u-bordeaux.fr>>>>>> wrote:
> >      >      >      >      >
> >      >      >      >      >     Hi all,
> >      >      >      >      >
> >      >      >      >      >     Can I make any assumption on the
> >     orientation of
> >      >     triangular
> >      >      >      >     facets in a
> >      >      >      >      >     tetrahedral plex ? I need the inward
> facet
> >      >     normals. Do
> >      >      >     I need
> >      >      >      >     to use
> >      >      >      >      >     DMPlexGetOrientedFace or can I rely on
> >     either
> >      >     the tet
> >      >      >     vertices
> >      >      >      >      >     ordering,
> >      >      >      >      >     or the faces ordering ? Could
> >      >      >     DMPlexGetRawFaces_Internal be
> >      >      >      >     enough ?
> >      >      >      >      >
> >      >      >      >      >
> >      >      >      >      > You can do it by hand, but you have to
> >     account for
> >      >     the face
> >      >      >      >     orientation
> >      >      >      >      > relative to the cell. That is what
> >      >      >      >      > DMPlexGetOrientedFace() does. I think it
> >     would be
> >      >     easier
> >      >      >     to use the
> >      >      >      >      > function below.
> >      >      >      >      >
> >      >      >      >      >     Alternatively, is there a function that
> >      >     computes the
> >      >      >     normals
> >      >      >      >     - without
> >      >      >      >      >     bringing out the big guns ?
> >      >      >      >      >
> >      >      >      >      >
> >      >      >      >      > This will compute the normals
> >      >      >      >      >
> >      >      >      >      >
> >      >      >      >
> >      >      >
> >      >
> >
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexComputeCellGeometryFVM.html
> >      >      >      >      > Should not be too heavy weight.
> >      >      >      >      >
> >      >      >      >      >    THanks,
> >      >      >      >      >
> >      >      >      >      >      Matt
> >      >      >      >      >
> >      >      >      >      >     Thanks
> >      >      >      >      >
> >      >      >      >      >     --
> >      >      >      >      >     Nicolas
> >      >      >      >      >
> >      >      >      >      >
> >      >      >      >      >
> >      >      >      >      > --
> >      >      >      >      > 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/~knepley/>
> >      >      >      >
> >      >      >      >
> >      >      >      >
> >      >      >      > --
> >      >      >      > 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/~knepley/>
> >      >      >
> >      >      >
> >      >      >
> >      >      > --
> >      >      > 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/~knepley/>
> >      >
> >      >
> >      >
> >      > --
> >      > 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/~knepley/>
> >
> >
> >
> > --
> > 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/~knepley/>
>


-- 
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/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210308/af3da83d/attachment-0001.html>


More information about the petsc-users mailing list