[petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM
Matthew Knepley
knepley at gmail.com
Mon Nov 14 17:00:01 CST 2022
On Mon, Nov 14, 2022 at 4:19 PM Blaise Bourdin <bourdin at mcmaster.ca> wrote:
> Replying to myself so that it gets into the googles:
>
> The reference tetrahedron used by DMPlexComputeCellGeometryAffineFEM has
> vertices at (-1,1,-1), (-1,-1,-1), (1, -1, -1) and (-1,-1,1).
>
It is
(-1, -1, -1) -- (-1, 1, -1) -- (1, -1, -1) -- (-1, -1, 1)
that way the first face has an outward normal.
Matt
> Blaise
>
>
> On Nov 10, 2022, at 6:42 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
> 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/>
>
>
> —
> 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/20221114/6699d237/attachment-0001.html>
More information about the petsc-users
mailing list