[petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM
Matthew Knepley
knepley at gmail.com
Thu Nov 10 17:42:48 CST 2022
On Thu, Nov 10, 2022 at 3:46 PM Blaise Bourdin <bourdin at mcmaster.ca> wrote:
> I am not sure I am buying this… If the tet was inverted, detJ would be
> negative, but it is always 1/8, as expected.
>
> The attached mesh is a perfectly valid tet generated by Cubit, with
> orientation matching the exodus documentation (ignore the mid-edge dof
> since this is a tet4).
> Here is what I get out of the code I attached in my previous email:
>
Yes, I use the opposite convention from ExodusII. In my opinion, orienting
face (1, 2, 3) to have an inward normal is sacrilegious.
Thanks,
Matt
> *SiMini*:Tests (dmplex)$ ./TestDMPlexComputeCellGeometryAffineFEM -i
> ../TestMeshes/TetCubit.gen
>
> filename ../TestMeshes/TetCubit.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
>
> From J, invJ, and v0, I still can’t reconstruct a reasonable reference tet
> which I was naively assuming was either the unit simplex, or the simplex
> with vertices at (-1,-1,-1), (-1,0,-1), (0, -1, -1), and (-1,-1,1) not
> necessarily in this order. In order to build my FE basis functions on the
> reference element, I really need to know what this element is…
>
> Blaise
>
>
>
>
> On Nov 9, 2022, at 6:56 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
> 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/>
>
>
> —
> 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/20221110/c057bd8a/attachment-0001.html>
More information about the petsc-users
mailing list