[petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Matthew Knepley]

On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin <bourdin at mcmaster.ca> wrote:

>
>
> On Nov 9, 2022, at 10:04 AM, Matthew Knepley <knepley at gmail.com> wrote:
>
> On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin <bourdin at mcmaster.ca> wrote:
>
> Hi,
>
> What reference simplex is DMPlexComputeCellGeometryAffineFEM using in 2
> and 3D?
> I am used to computing my shape functions on the unit simplex (vertices at
> the origin and each e_i), but it does not look to be the reference simplex
> in this function:
>
> In 3D, for the unit simplex with vertices at (0,0,0) (1,0,0) (0,1,0)
> (0,0,1) (in this order), I get J = 1 / 2 . [[-1,-1,-1],[1,0,0],[0,0,1]] and
> v0 = [0,0,1]
>
> In 2D, for the unit simplex with vertices at (0,0), (1,0), and (0,1), I
> get J = 1 / 2. I and v0 = [0,0], which does not make any sense to me (I was
> assuming that the 2D reference simplex had vertices at (-1,-1), (1, -1) and
> (-1,1), but if this were the case, v0 would not be 0).
>
> I can build a simple example with meshes consisting only of the unit
> simplex in 2D and 3D if that would help.
>
>
> I need to rewrite the documentation on geometry, but I was waiting until I
> rewrite the geometry calculations to fit into libCEED. Toby found a nice
> way to express them in BLAS form which I need to push through everything.
>
> I always think of operating on the cell with the first vertex at the
> origin (I think it is easier), so I have a xi0 that translates the first
> vertex
> of the reference to the origin, and a v0 that translates the first vertex
> of the real cell to the origin. You can see this here
>
>
> https://gitlab.com/petsc/petsc/-/blob/main/include/petsc/private/petscfeimpl.h#L251
>
> This explains the 2D result. I cannot understand your 3D result, unless
> the vertices are in another order.
>
>
> That makes two of us, then… I am attaching a small example and test meshes
> (one cell being the unit simplex starting with the origin and numbered in
> direct order when looking from (1,1,1)
>

Oh, it is probably inverted. All faces are oriented for outward normals. It
is in the Orientation chapter in the book :)

  Thanks,

      Matt


>  filename ../TestMeshes/1Tri.gen
>
> Vec Object: coordinates 1 MPI process
>
>   type: seq
>
> 0.
>
> 0.
>
> 1.
>
> 0.
>
> 0.
>
> 1.
>
> v0
>
>  0:   0.0000e+00   0.0000e+00
>
> J
>
>  0:   5.0000e-01   0.0000e+00
>
>  0:   0.0000e+00   5.0000e-01
>
> invJ
>
>  0:   2.0000e+00  -0.0000e+00
>
>  0:  -0.0000e+00   2.0000e+00
>
> detJ : 0.25
>
> And
>
>  filename ../TestMeshes/1Tet.gen
>
> Vec Object: coordinates 1 MPI process
>
>   type: seq
>
> 0.
>
> 0.
>
> 0.
>
> 1.
>
> 0.
>
> 0.
>
> 0.
>
> 1.
>
> 0.
>
> 0.
>
> 0.
>
> 1.
>
> v0
>
>  0:   1.0000e+00   0.0000e+00   0.0000e+00
>
> J
>
>  0:  -5.0000e-01  -5.0000e-01  -5.0000e-01
>
>  0:   5.0000e-01   0.0000e+00   0.0000e+00
>
>  0:   0.0000e+00   0.0000e+00   5.0000e-01
>
> invJ
>
>  0:   0.0000e+00   2.0000e+00   0.0000e+00
>
>  0:  -2.0000e+00  -2.0000e+00  -2.0000e+00
>
>  0:   0.0000e+00   0.0000e+00   2.0000e+00
>
> detJ : 0.125
>
> I don’t understand why v0=(0,0) in 2D and (1,0,0) in 3D (but don’t really
> care) since I only want J. J makes no sense to me in 3D. In particular, one
> does not seem to have X~ = invJ.X + v0 (X = J.(X~-v0) as stated in
> CoordinatesRefToReal (it works in 2D if V0 = (1,1), which is consistent
> with a reference simplex with vertices at (-1,-1), (1,-1) and (-1,1)).
>
> What am I missing?
>
> Blaise
>
> /
>
>
>
>   Thanks,
>
>       Matt
>
>
> Regards,
> Blaise
>
>
>
> —
> Canada Research Chair in Mathematical and Computational Aspects of Solid
> Mechanics (Tier 1)
> Professor, Department of Mathematics & Statistics
> Hamilton Hall room 409A, McMaster University
> 1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada
> https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243
>
>
>
> --
> 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/>
>
>
> —
> Canada Research Chair in Mathematical and Computational Aspects of Solid
> Mechanics (Tier 1)
> Professor, Department of Mathematics & Statistics
> Hamilton Hall room 409A, McMaster University
> 1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada
> https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243
>
>

-- 
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/20221109/57c35e72/attachment.html>