Pyramids #623

lijas · 2023-03-20T10:02:27Z

Adds RefPyramid, Pyramid cells, Lagrange{RefPyramid, 1&2}

fredrikekre · 2023-03-20T11:18:53Z

src/interpolations.jl

+    zzero = z ≈ 1.0
+    i == 1 && return zzero ? 0.0 : (-x*y+(z-1)*(-x-y-z+1))/(z-1)


This means the derivative can't be evaluated at z = 1.0 (truncating the dual part of the input, so it will have 0 derivative).

At least not until ForwardDiff 0.16.

But this is discontinuous around z = 1? Is that intended?

I used this as ref https://defelement.com/elements/examples/pyramid-Lagrange-1.html
I think that the pyramid is degenerate at the fifth node. But note that N_I = 1 at node I, and zero at the other nodes is still true for pyramids

Maybe @termi-official has some insight regarding pyramid elements

One of the basis functions is $-xy/(z-1)$. So yes, you have a singularity at $z=1$ (note $x=y=0$). IIRC this is a problem with the way how these pyramidal elements are constructed. ~~If you want the derivative at this specific point you can formulate y and x as a function of z and apply l'Hospital.~~

Wieners, C. (1997). Conforming discretizations on tetrahedrons, pyramids, prisms and hexahedrons. Preprint, University of Stuttgart. has a proof that there exists no continuous differentiable (bi-)linear function on pyramids.

fredrikekre · 2023-03-21T12:12:23Z

src/interpolations.jl

+    i == 2 && return zzero ? 0.0 : x*(y+z-1)/(z-1)
+    i == 3 && return zzero ? 0.0 : -x*y/(z-1)
+    i == 4 && return zzero ? 0.0 : y*(x+z-1)/(z-1)
+    i == 5 && return zzero ? 1.0 : z


This one is at least not problematic:

Suggested change

i == 5 && return zzero ? 1.0 : z

i == 5 && return z

For the other ones, if we are sufficiently close to z = 1, can we replac them with expressions that have the same value and derivative in the limit perhaps?

There is no derivative at and close to z=1 with x=y=0. See e.g. l'Hospital on the basis function i=3 with the coordinates x and y depending on z gives you $-(x'*y + y'*x)/(z'-1) = 0/0$ again. Maybe we should throw an error in this case for the derivative?

termi-official · 2023-06-30T07:17:48Z

src/docs.jl

+Vertex numbers:                                  | Vertex coordinates:
+                                                 | v1: 𝛏 = (-1.0, -1.0, -1.0)
+         5 -. _               5 -. _             | v2: 𝛏 = ( 1.0, -1.0, -1.0)
+         | \    \._           |\     \._         | v3: 𝛏 = ( 1.0,  1.0, -1.0)
+  ^ ξ₃   |   \      \ _       | \        \ _     | v4: 𝛏 = (-1.0,  1.0, -1.0)
+  |  ξ₂  |     \       3      |  4--------- 3    | v5: 𝛏 = (-1.0, -1.0,  1.0)
+  | /    |       \    /       | /          /     | v6: 𝛏 = ( 1.0, -1.0,  1.0)
+  |/_ _  |         \ /        |/          /      | v7: 𝛏 = ( 1.0,  1.0,  1.0)
+      ξ₁ 1----------2         1----------2       | v8: 𝛏 = (-1.0,  1.0,  1.0)


I think the coordinates are wrong.

termi-official · 2023-06-30T07:19:01Z

src/docs.jl

+            +----8---+        +----8---+ |
+          5/        /|      5/|        | |  e1: (v1, v2),  e2: (v2, v3)
+          /       7/ |12    / |9     12| |  e3: (v3, v4),  e4: (v4, v1)
+         +----6---+  |     +  |        | |  e5: (v5, v6),  e6: (v6, v7)
+         |        |  +     |  +---4----+ |  e7: (v7, v8),  e8: (v8, v5)
+       10|      11| /    10| /1       /  |  e9: (v1, v5), e10: (v2, v6)
+         |        |/3      |/        /3  | e11: (v3, v7), e12: (v4, v8)
+         +---2----+        +---2----+    |


I think the edge numbering and picture is wrong.

lijas · 2023-06-30T07:50:10Z

I have checked the heat_equation.jl with pyramid elements and it converges to the analytical solution (notably a bit slower than with for example Hexahedrons). It seems like PR Helper script for convergence rates for Poisson problems #640 will check this
PR More reference cell documentation #689 will update the docs for all reference shapes/cells, so perhaps there is no need to add docs about the RefPyramid in this PR?
Regarding the singularity at the 5th node. This does not seems to be an issue in practice, because the derivative of the shape functions never evaluated at the nodes, right? Perhaps if you have some :lobatto integration scheme...

lijas · 2023-06-30T09:40:23Z

I added quadratic interpolations aswell. I also noticed that I had not used defelements numberings scheme, so I updated that aswell. (Note however that did not match VTK numbering scheme for pyramids)

The quadratic elements does not coverge for the heat equations problem, so there is something wrong...

codecov-commenter · 2023-06-30T09:47:03Z

Codecov Report

Patch coverage: 73.18% and project coverage change: -0.53% ⚠️

Comparison is base (fd4cfd3) 92.40% compared to head (c94d19b) 91.88%.

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #623      +/-   ##
==========================================
- Coverage   92.40%   91.88%   -0.53%     
==========================================
  Files          32       33       +1     
  Lines        4779     4916     +137     
==========================================
+ Hits         4416     4517     +101     
- Misses        363      399      +36

Files Changed	Coverage Δ
src/Ferrite.jl	`100.00% <ø> (ø)`
src/Quadrature/gaussquad_pyramid_table.jl	`0.00% <0.00%> (ø)`
src/Quadrature/quadrature.jl	`85.43% <8.33%> (-10.22%)`	⬇️
src/Export/VTK.jl	`93.54% <100.00%> (+1.88%)`	⬆️
src/Grid/grid.jl	`92.43% <100.00%> (+0.34%)`	⬆️
src/Grid/grid_generators.jl	`100.00% <100.00%> (ø)`
src/interpolations.jl	`97.01% <100.00%> (+0.26%)`	⬆️

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

lijas · 2023-07-03T15:08:56Z

The quadratic interpolation for pyramids does not converge for some reason, even though in passes all the tests in test_interpolation.jl. I have double checked with defelements, and I am using the same numbering scheme.

@termi-official Any idea what the issue could be? 🤔

termi-official · 2023-07-03T19:10:24Z

The quadrature rule seems to be broken. Choose N=20 and qr order=3 for the linear element and it also diverges.

lijas · 2023-07-04T09:28:43Z

Ok thanks for checking. Is N=20 the number of elements in each direction? Maybe qr order=3 is simply to few points?
Currently the pyramids are using the same type of quadrature rule as the wedges (:polyquad). Here is another paper deriving quadrature rule for pyramids. I will try to use these points instead to see if it works better.

Edit: From the paper: "We note that pyramidal elements are unique among the family of three-dimensional finite elements in that the functions spaces used span complete polynomials of a fixed degree and also include nonpolynomial functions; see, for example, [14]. Thus, even for relatively low-order pyramidal elements, sufficiently high-degree integration formulas should be used to minimize approximation error."

termi-official · 2023-07-04T12:27:24Z

Is N=20 the number of elements in each direction?

Yes.

Maybe qr order=3 is simply to few points?

As stated in your edit, there is no exact quadrature rule for pyramids, because some basis functions are fractional polynomials. This was just an example. If you take the quadratic pyramids and increase the quadrature, then the error should not approximately double, even if you are not exactly integrating. The solution to this heat problem is actually super simple and we use a regular discretization, so there is no reason that any of your basis functions should heavily oscillate on such a fine mesh on any element.

Currently the pyramids are using the same type of quadrature rule as the wedges (:polyquad). Here is another paper deriving quadrature rule for pyramids. I will try to use these points instead to see if it works better.

I do not think that polyquad is per-se the issue, because it is known to work well. However, I do not fully understand your transformation in the quadrature rule. You seem to move some of the integration points outside the integration domain, which breaks the numerical integration.

julia> qr = QuadratureRule{RefPyramid}(2)
QuadratureRule{RefPyramid, Float64, 3}([0.6028728035309393, 0.5159484657839318, 0.5159484657839318, 0.5159484657839318, 0.5159484657839318], Vec{3, Float64}[[0.39170896144022116, 0.39170896144022116, 0.6082910385597788], [1.5735845722389405, 0.8546635164271443, 0.14533648357285572], [0.8546635164271443, 1.5735845722389405, 0.14533648357285572], [0.13574246061534806, 0.8546635164271443, 0.14533648357285572], [0.8546635164271443, 0.13574246061534806, 0.14533648357285572]])

julia> qr.points
5-element Vector{Vec{3, Float64}}:
 [0.39170896144022116, 0.39170896144022116, 0.6082910385597788]
 [1.5735845722389405, 0.8546635164271443, 0.14533648357285572]
 [0.8546635164271443, 1.5735845722389405, 0.14533648357285572]
 [0.13574246061534806, 0.8546635164271443, 0.14533648357285572]
 [0.8546635164271443, 0.13574246061534806, 0.14533648357285572]

lijas · 2023-07-04T14:37:47Z

Thanks @termi-official . There was in issue with the mapping of the quadrature points. Now it seems like the both the linear and quadratic interpolations converge.

Do you have any more feedback on this PR? Do you think we should add some docs in this PR or #689 ?

termi-official

Do you have any more feedback on this PR? Do you think we should add some docs in this PR or #689 ?

Since the linked PR will take care of the docs (and seems to have already for Pyramids) lets leave it separate.

KnutAM · 2023-07-04T15:15:34Z

src/interpolations.jl

+function shape_value(ip::Lagrange{RefPyramid,1}, ξ::Vec{3}, i::Int)
+    (x,y,z) = ξ
+    zzero = z ≈ 1.0
+    i == 1 && return zzero ? 0.0 : (-x*y+(z-1)*(-x-y-z+1))/(z-1)


Suggested change

i == 1 && return zzero ? 0.0 : (-x*y+(z-1)*(-x-y-z+1))/(z-1)

i == 1 && return zzero ? zero(x) : (-x*y+(z-1)*(-x-y-z+1))/(z-1)

To make type-stable for non-Float64 input? (Same below for 2nd order)

lijas · 2023-07-27T12:57:11Z

Docs will be fixed in another PR, but in case we want ascii-documentation in the future, I paste possible representations here:

alt 1

Vertex numbers:                                  | Vertex coordinates:
                                                 |
         5 -. _               5 -. _             |
         | \    \._           |\     \._         | 
  ^ ξ₃   |   \      \ _       | \        \ _     | 
  |  ξ₂  |     \       3      |  4--------- 3    | 
  | /    |       \    /       | /          /     | 
  |/_ _  |         \ /        |/          /      |
      ξ₁ 1----------2         1----------2       |

Alt 2:

Vertex numbers:                                  | Vertex coordinates:
         5 _                 5 _           
        / \  \ _            /|   \ _         
       .   |     \         . |       \        
^ ξ₃    |    |      4      | 1--------4      
|      .       |   /      . /        /        
+-> ξ₂ |        \ /       |/        /
/      2---------3        2--------3           
ξ₁

fredrikekre reviewed Mar 20, 2023

View reviewed changes

fredrikekre reviewed Mar 21, 2023

View reviewed changes

fredrikekre added this to the 0.4.0 milestone Mar 23, 2023

termi-official linked an issue Mar 23, 2023 that may be closed by this pull request

Support for Wedge&Pyramid elements #580

Closed

2 tasks

termi-official mentioned this pull request Mar 24, 2023

Support for Wedge&Pyramid elements #580

Closed

2 tasks

lijas added 2 commits June 29, 2023 14:17

fake rebase to master

704a1eb

wip tests and docs

07a9473

lijas force-pushed the eb_pyramid2 branch from 8ded277 to 07a9473 Compare June 29, 2023 13:39

lijas added 2 commits June 29, 2023 17:51

add more tests for pyramid

01cb8dc

add checksum for wedge and pyramid

a63be10

termi-official reviewed Jun 30, 2023

View reviewed changes

clean up some files

3363eca

add quadratic interpolation and fix numberings w.r.t defelemnts

b5d9312

lijas added 3 commits June 30, 2023 11:52

add more gauss points for pyramids

bbf5435

fix face and edge numbering for pyramid

2caecd0

clean up

3bf4db8

fix mapping of quadrature rule

61ab4a5

termi-official approved these changes Jul 4, 2023

View reviewed changes

KnutAM reviewed Jul 4, 2023

View reviewed changes

lijas added 2 commits July 27, 2023 13:44

use abstract float type

632029f

Merge branch 'master' into eb_pyramid2

c94d19b

fredrikekre merged commit d474caf into Ferrite-FEM:master Jul 27, 2023

fredrikekre added the needs changelog label Jul 27, 2023

KnutAM mentioned this pull request Jul 30, 2023

Support wedge and pyramid elements Ferrite-FEM/FerriteMeshParser.jl#21

Open

		zzero = z ≈ 1.0
		i == 1 && return zzero ? 0.0 : (-xy+(z-1)(-x-y-z+1))/(z-1)

	i == 1 && return zzero ? 0.0 : (-xy+(z-1)(-x-y-z+1))/(z-1)
	i == 1 && return zzero ? zero(x) : (-xy+(z-1)(-x-y-z+1))/(z-1)

Pyramids #623

Pyramids #623

Uh oh!

Conversation

lijas commented Mar 20, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

termi-official Mar 20, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lijas commented Jun 30, 2023

Uh oh!

lijas commented Jun 30, 2023

Uh oh!

codecov-commenter commented Jun 30, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

lijas commented Jul 3, 2023

Uh oh!

termi-official commented Jul 3, 2023

Uh oh!

lijas commented Jul 4, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

termi-official commented Jul 4, 2023

Uh oh!

lijas commented Jul 4, 2023

Uh oh!

termi-official left a comment

Choose a reason for hiding this comment

Uh oh!

KnutAM Jul 4, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lijas commented Jul 27, 2023

Uh oh!

Uh oh!

lijas commented Mar 20, 2023 •

edited

Loading

termi-official Mar 20, 2023 •

edited

Loading

codecov-commenter commented Jun 30, 2023 •

edited

Loading

lijas commented Jul 4, 2023 •

edited

Loading

KnutAM Jul 4, 2023 •

edited

Loading