-
Notifications
You must be signed in to change notification settings - Fork 90
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) #790
base: master
Are you sure you want to change the base?
Conversation
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 ReportPatch coverage:
❗ 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
☔ View full report in Codecov by Sentry. |
…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
return (face1, face2, face3, face4) | ||
for (name, facefuncs) in ( | ||
(:ArbitraryOrderLagrange, :facedof_indices), | ||
(:ArbitraryOrderDiscontinuousLagrange, :dirichlet_facedof_indices) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@fredrikekre should we separate the fixed order implementation which we have from the arbitrary order infrastructure?
src/interpolations.jl
Outdated
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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. :)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Passing Val(order)
might be a quick and easy remedy
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
# Based on vtkTriQuadraticHexahedron (see https://kitware.github.io/vtk-examples/site/Cxx/GeometricObjects/IsoparametricCellsDemo/) | ||
# Permutation to switch numbering to Ferrite ordering |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can you elaborate?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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)
src/interpolations.jl
Outdated
@@ -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}` |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we should wrap this also into the arbitrary order Lagrange struct.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
Just as a reminder, this should come after #829 right? |
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 |
For #626
Some comments:
Note: I also edited the existing Tri345 to be arbitrary but it doesn't have the arbitrary basis/points thing