Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Matthew Knepley]

On Mon, Nov 14, 2022 at 4:19 PM Blaise Bourdin  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  wrote:
>
> On Thu, Nov 10, 2022 at 3:46 PM Blaise Bourdin 
> 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.e+00   0.e+00   0.e+00
>> J
>>  0:  -5.e-01  -5.e-01  -5.e-01
>>  0:   5.e-01   0.e+00   0.e+00
>>  0:   0.e+00   0.e+00   5.e-01
>> invJ
>>  0:   0.e+00   2.e+00   0.e+00
>>  0:  -2.e+00  -2.e+00  -2.e+00
>>  0:   0.e+00   0.e+00   2.e+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  wrote:
>>
>> On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin 
>> wrote:
>>
>>
>>
>> On Nov 9, 2022, at 10:04 AM, Matthew Knepley  wrote:
>>
>> On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin 
>> 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.e+00   0.e+00
>> J
>>  0:   5.e-01   0.e+00
>>  0:   0.e+00   5.e-01
>> invJ
>>  0:   2.e+00  -0.e+00
>>  0:  -0.e+00   2.e+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.e+00   0.e+00   0.e+00
>> J
>>  0:  -5.e-01  -5.e-01  -5.e-01
>>  0:   5.e-01   0.e+00   0.e+00
>>  0:   0.e+00   0.e+00   5.e-01

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Blaise Bourdin]

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  wrote:



On Thu, Nov 10, 2022 at 3:46 PM Blaise Bourdin  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.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+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  wrote:



On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin  wrote:







On Nov 9, 2022, at 10:04 AM, Matthew Knepley  wrote:



On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin  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.e+00   0.e+00

J

 0:   5.e-01   0.e+00

 0:   0.e+00   5.e-01

invJ

 0:   2.e+00  -0.e+00

 0:  -0.e+00   2.e+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.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Matthew Knepley]

On Thu, Nov 10, 2022 at 3:46 PM Blaise Bourdin  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.e+00   0.e+00   0.e+00
>
> J
>
>  0:  -5.e-01  -5.e-01  -5.e-01
>
>  0:   5.e-01   0.e+00   0.e+00
>
>  0:   0.e+00   0.e+00   5.e-01
>
> invJ
>
>  0:   0.e+00   2.e+00   0.e+00
>
>  0:  -2.e+00  -2.e+00  -2.e+00
>
>  0:   0.e+00   0.e+00   2.e+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  wrote:
>
> On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin 
> wrote:
>
>
>
> On Nov 9, 2022, at 10:04 AM, Matthew Knepley  wrote:
>
> On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin  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.e+00   0.e+00
> J
>  0:   5.e-01   0.e+00
>  0:   0.e+00   5.e-01
> invJ
>  0:   2.e+00  -0.e+00
>  0:  -0.e+00   2.e+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.e+00   0.e+00   0.e+00
> J
>  0:  -5.e-01  -5.e-01  -5.e-01
>  0:   5.e-01   0.e+00   0.e+00
>  0:   0.e+00   0.e+00   5.e-01
> invJ
>  0:   0.e+00   2.e+00   0.e+00
>  0:  -2.e+00  -2.e+00  -2.e+00
>  0:   0.e+00   0.e+00   2.e+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
>
> /
>

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Blaise Bourdin]

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:




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.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+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  wrote:



On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin  wrote:







On Nov 9, 2022, at 10:04 AM, Matthew Knepley  wrote:



On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin  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.e+00   0.e+00

J

 0:   5.e-01   0.e+00

 0:   0.e+00   5.e-01

invJ

 0:   2.e+00  -0.e+00

 0:  -0.e+00   2.e+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.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+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 

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Matthew Knepley]

On Wed, Nov 9, 2022 at 10:46 AM Blaise Bourdin  wrote:

>
>
> On Nov 9, 2022, at 10:04 AM, Matthew Knepley  wrote:
>
> On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin  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.e+00   0.e+00
>
> J
>
>  0:   5.e-01   0.e+00
>
>  0:   0.e+00   5.e-01
>
> invJ
>
>  0:   2.e+00  -0.e+00
>
>  0:  -0.e+00   2.e+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.e+00   0.e+00   0.e+00
>
> J
>
>  0:  -5.e-01  -5.e-01  -5.e-01
>
>  0:   5.e-01   0.e+00   0.e+00
>
>  0:   0.e+00   0.e+00   5.e-01
>
> invJ
>
>  0:   0.e+00   2.e+00   0.e+00
>
>  0:  -2.e+00  -2.e+00  -2.e+00
>
>  0:   0.e+00   0.e+00   2.e+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/
> 
>
>
> —
> 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/ 

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Blaise Bourdin]

On Nov 9, 2022, at 10:04 AM, Matthew Knepley  wrote:



On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin  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)


 filename ../TestMeshes/1Tri.gen

Vec Object: coordinates 1 MPI process

  type: seq

0.

0.

1.

0.

0.

1.

v0

 0:   0.e+00   0.e+00

J

 0:   5.e-01   0.e+00

 0:   0.e+00   5.e-01

invJ

 0:   2.e+00  -0.e+00

 0:  -0.e+00   2.e+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.e+00   0.e+00   0.e+00

J

 0:  -5.e-01  -5.e-01  -5.e-01

 0:   5.e-01   0.e+00   0.e+00

 0:   0.e+00   0.e+00   5.e-01

invJ

 0:   0.e+00   2.e+00   0.e+00

 0:  -2.e+00  -2.e+00  -2.e+00

 0:   0.e+00   0.e+00   2.e+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/



























— 
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






















TestDMPlexComputeCellGeometryAffineFEM.c
Description: TestDMPlexComputeCellGeometryAffineFEM.c


1Tet.gen
Description: 1Tet.gen


1Tri.gen
Description: 1Tri.gen

Re: [petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Matthew Knepley]

On Tue, Nov 8, 2022 at 9:14 PM Blaise Bourdin  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.

  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/ 

[petsc-users] Reference element in DMPlexComputeCellGeometryAffineFEM

[Blaise Bourdin]

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.

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