On Mon, Mar 8, 2021 at 11:18 AM Nicolas Barral < [email protected]> wrote:
> On 08/03/2021 15:55, Matthew Knepley wrote: > > On Mon, Mar 8, 2021 at 4:02 AM Nicolas Barral > > <[email protected] > > <mailto:[email protected]>> wrote: > > > > On 07/03/2021 22:56, Matthew Knepley wrote: > > > On Sun, Mar 7, 2021 at 4:51 PM Nicolas Barral > > > <[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>> wrote: > > > > > > > > > On 07/03/2021 22:30, Matthew Knepley wrote: > > > > On Sun, Mar 7, 2021 at 4:13 PM Nicolas Barral > > > > <[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>> > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>>> wrote: > > > > > > > > On 07/03/2021 16:54, Matthew Knepley wrote: > > > > > On Sun, Mar 7, 2021 at 8:52 AM Nicolas Barral > > > > > <[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>> > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>> > > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>> > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>>>> 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 > > > > > > <[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>> > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>> > > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>> > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>>> > > > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>> > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>> > > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>> > > > > <mailto:[email protected] > > <mailto:[email protected]> > > > <mailto:[email protected] > > <mailto:[email protected]>>>>>> 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/>
