Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) #790

AbdAlazezAhmed · 2023-08-14T22:54:28Z

For #626
Some comments:

using Vandermonde inverse resulted in worse floating point errors during Dirac-delta tests
Tested order 10 with heat equation example, seems to work (haven't checked convergence tho but seems ok from paraview)

Note: I also edited the existing Tri345 to be arbitrary but it doesn't have the arbitrary basis/points thing

termi-official · 2023-08-14T22:58:37Z

We definitely should be able to pass 1D basis functions at least for hypercubes, with GLL as default for continuous Lagrange and GL as default for discontinuous Lagrange (for the high order infrastructure), because we can rule out Runge's phenomenon.

why? Best regards.

codecov-commenter · 2023-08-15T14:57:38Z

Codecov Report

Patch coverage: 98.82% and project coverage change: +0.40% 🎉

Comparison is base (b7f15e4) 91.88% compared to head (17040e1) 92.28%.
Report is 2 commits behind head on master.

❗ 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     #790      +/-   ##
==========================================
+ Coverage   91.88%   92.28%   +0.40%     
==========================================
  Files          33       33              
  Lines        4916     5198     +282     
==========================================
+ Hits         4517     4797     +280     
- Misses        399      401       +2

Files Changed	Coverage Δ
src/interpolations.jl	`97.64% <98.82%> (+0.63%)`	⬆️

... and 1 file with indirect coverage changes

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

…move them and use the already implemented ones. Performance isn't too bad, values don't allocate since everything is in the constructor

order 13 hex tests from 1m39.7 sec to 1m35.3sec

AbdAlazezAhmed · 2023-08-19T11:59:04Z

src/interpolations.jl

-    return (face1, face2, face3, face4)
+for (name,  facefuncs) in (
+    (:ArbitraryOrderLagrange,               :facedof_indices),
+    (:ArbitraryOrderDiscontinuousLagrange,  :dirichlet_facedof_indices)


I just realized that using GL points should result in dirichlet_facedof_indices being empty IIUC. This makes this function need to either check for faces on boundary or have predetermined options for the points (equidistant, GL, Radau) instead of passing the points in the constructor.
I think having predetermined options would be easier and more consistent (as we guarantee the ordering), what do you think?

I would suggest that we just throw an error for now when trying to extract the Dirichlet dofs for all L2 elements which are not based on "symmetric" nodes, when the basis functions of the interior nodes are non-zero on the boundary (and for modal basis functions). If we want to work further on strong enforcement of the Dirichlet condition we have to think a bit how to do it correctly.

termi-official · 2023-08-19T21:50:03Z

src/interpolations.jl

@fredrikekre should we separate the fixed order implementation which we have from the arbitrary order infrastructure?

termi-official · 2023-08-19T21:55:29Z

src/interpolations.jl

+function getPermLagrangeTri(order)
+    verts = SVector{3}(order + 1, (order + 1) * (order + 2) ÷ 2, 1)
+    edge1 = SVector{order-1}((sum(i:order+1) for i in order:-1:2))
+    edge2 = SVector{order-1}((sum(i:order+1)+1 for i in 3:(order+1)))
+    edge3 = SVector{order-1}((i for i in 2:order))
+    # TODO: fix this bottleneck, nested generators don't work so one has to use ...
+    interior = NTuple{(order - 2) * (order - 1) ÷ 2}((i for j in (order+1):-1:4 for i in sum(j:order+1)+2 : sum(j:order+1)+j-2))
+    return SVector((verts..., edge1..., edge2..., edge3..., interior...))
+end


I think it is worth to investigate a code generator here in the future (e.g. to just generate the elements from the Ciarlet definition). So, no need to invest too much time here if it is not absurdely slow. :)

Passing Val(order) might be a quick and easy remedy

termi-official

Comments below. I think we need to add some doc strings and think about discoverability. Not sure if we can study the convergence order at this point.

termi-official · 2023-08-19T22:07:07Z

src/interpolations.jl

+for (name,  default_coords) in (
+    (:ArbitraryOrderLagrange,               :both),
+    (:ArbitraryOrderDiscontinuousLagrange,  :neither)
+    )
+    @eval begin
+        function $(name){RefQuadrilateral, order}(points::Vector{Float64} = GaussQuadrature.legendre(order+1, GaussQuadrature.$(default_coords))[1]) where order


Can we add doc strings for these, so users know about the functionality?

Also, how can users discover this feature? Maybe we should add some hint in the Lagrange doc string and a tip one of the examples? Not sure if it is worth to copy paste the full heat example here. In the future we can think about adding a Lagrangian formulation of some advection problem.

termi-official · 2023-08-19T22:08:34Z

src/interpolations.jl

+# Based on vtkTriQuadraticHexahedron (see https://kitware.github.io/vtk-examples/site/Cxx/GeometricObjects/IsoparametricCellsDemo/)
+# Permutation to switch numbering to Ferrite ordering


Can you elaborate?

The first line refers to the numbering convention used (copied from fixed order one)
The second one refers to the permutation that maps from the generated (-1:1)^n coordinates to the ones following the convention (such as {-1,...,1} ->{-1,1,...} in 1D)

termi-official · 2023-08-26T09:32:08Z

src/interpolations.jl

@@ -32,6 +29,17 @@ The following interpolations are implemented:
 * `Serendipity{RefQuadrilateral,2}`
 * `Serendipity{RefHexahedron,2}`

+The following interpolations are implemented with arbitrary order:
+
+* `Lagrange{RefTriangle,order}`


I think we should wrap this also into the arbitrary order Lagrange struct.

I think the name is kinda too long 😄. Maybe we can make it a fallback for Lagrange constructor so we don't have to use this long name?

termi-official · 2023-11-15T14:57:15Z

Just as a reminder, this should come after #829 right?

AbdAlazezAhmed · 2023-11-15T15:19:05Z

Just as a reminder, this should come after #829 right?

Yeah I think so as #829 would eliminate the need for permutations IIUC.

termi-official · 2024-01-04T23:53:45Z

I think this one here might be helpful for future reference on simplices: Warburton, T. An explicit construction of interpolation nodes on the simplex. J Eng Math 56, 247–262 (2006). https://doi.org/10.1007/s10665-006-9086-6

AbdAlazezAhmed added 2 commits August 14, 2023 18:54

Initial commit, bad code

681d93b

Draft

aea12d0

Dear compiler,

84f9dc0

why? Best regards.

AbdAlazezAhmed added 8 commits August 17, 2023 02:53

Line + Quad + start of Hex + some GL/GLL generators, will probably re…

cdcf12f

…move them and use the already implemented ones. Performance isn't too bad, values don't allocate since everything is in the constructor

Now hex works with bad performance

9d8fc66

a bit faster but is it worth it? (less allocations)

324a732

order 13 hex tests from 1m39.7 sec to 1m35.3sec

Don't reinvent the wheel :D

ff4347b

use existing GLL/GL code + some tests

b41a072

Discontinuous + some tests

f76cb72

another forgotten method

a21469a

some thing I forgot while copy-pasting

e9cbc21

AbdAlazezAhmed changed the title ~~[WIP] : Arbitrary order Lagrange interpolations~~ Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) Aug 18, 2023

AbdAlazezAhmed marked this pull request as ready for review August 18, 2023 21:56

AbdAlazezAhmed added enhancement gsoc23 Google Summer of Code 2023 labels Aug 18, 2023

AbdAlazezAhmed commented Aug 19, 2023

View reviewed changes

termi-official reviewed Aug 19, 2023

View reviewed changes

AbdAlazezAhmed added 7 commits August 22, 2023 14:53

better performance for higher order Lagrange trianlge

6ab1dbf

add errors

15af785

Some docs

0cf8d93

errors + tests

0ea1c67

typo

5e6818d

add hint in heat equation example

08163a4

1.6 test thing

b9c14c9

termi-official reviewed Aug 26, 2023

View reviewed changes

Lagrange{RefTriangle} -> ArbitraryOrderLagrange{RefTriangle}

8c24f7a

update heat equation tip

17040e1

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

Generic arbitrary order interpolations #626

Open

termi-official mentioned this pull request Oct 12, 2023

Batch evaluations #819

Merged

2 tasks

AbdAlazezAhmed mentioned this pull request Oct 18, 2023

Make DofHandler work with other(any?) dof distribution #829

Draft

AbdAlazezAhmed marked this pull request as draft November 15, 2023 15:40

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) #790

Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) #790

AbdAlazezAhmed commented Aug 14, 2023 •

edited

termi-official commented Aug 14, 2023

codecov-commenter commented Aug 15, 2023 •

edited

AbdAlazezAhmed Aug 19, 2023

termi-official Aug 19, 2023

termi-official Aug 19, 2023

termi-official Aug 19, 2023

KnutAM Aug 20, 2023 •

edited

termi-official left a comment

termi-official Aug 19, 2023

termi-official Aug 19, 2023

AbdAlazezAhmed Aug 20, 2023

termi-official Aug 26, 2023

AbdAlazezAhmed Aug 27, 2023

termi-official commented Nov 15, 2023

AbdAlazezAhmed commented Nov 15, 2023

termi-official commented Jan 4, 2024

		# Based on vtkTriQuadraticHexahedron (see https://kitware.github.io/vtk-examples/site/Cxx/GeometricObjects/IsoparametricCellsDemo/)
		# Permutation to switch numbering to Ferrite ordering

Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) #790

Are you sure you want to change the base?

Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) #790

Conversation

AbdAlazezAhmed commented Aug 14, 2023 • edited

termi-official commented Aug 14, 2023

codecov-commenter commented Aug 15, 2023 • edited

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KnutAM Aug 20, 2023 • edited

Choose a reason for hiding this comment

termi-official left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

termi-official commented Nov 15, 2023

AbdAlazezAhmed commented Nov 15, 2023

termi-official commented Jan 4, 2024

AbdAlazezAhmed commented Aug 14, 2023 •

edited

codecov-commenter commented Aug 15, 2023 •

edited

KnutAM Aug 20, 2023 •

edited