|
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).
Blaise
On Nov 10, 2022, at 6:42 PM, Matthew Knepley <knep...@gmail.com> 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
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
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
—
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
—
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
—
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
|
