[petsc-users] DMPlex tetrahedra facets orientation
Nicolas Barral
nicolas.barral at math.u-bordeaux.fr
Mon Mar 8 10:18:22 CST 2021
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 ?
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/>
More information about the petsc-users
mailing list