On Thu, Nov 10, 2022 at 3:46 PM Blaise Bourdin <bour...@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 <knep...@gmail.com> wrote: > > On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin <bour...@mcmaster.ca> > wrote: > > > > On Nov 9, 2022, at 10:04 AM, Matthew Knepley <knep...@gmail.com> wrote: > > On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin <bour...@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/>