WIP extend Lagrange elements #482
Conversation
|
Thanks for the PR! I think you have just missed the "cell dofs", i.e. dofs that are internal to the cell. (Edit: In practice this meant that you only assembled the first 9 rows/colums of the local matrix even though the full 10x10 matrix was computed correctly. This should have errored IMO, so I opened a PR to check that the size of the matrix matches the number of dofs: #483) Your test function passes with diff --git a/src/extend_interpolations_Ferriteorder.jl b/src/extend_interpolations_Ferriteorder.jl
index 4682c693..781c60db 100644
--- a/src/extend_interpolations_Ferriteorder.jl
+++ b/src/extend_interpolations_Ferriteorder.jl
@@ -3,8 +3,9 @@ function Ferrite.getnbasefunctions(::Lagrange{2,RefTetrahedron,order}) where ord
end
Ferrite.nvertexdofs(::Lagrange{2,RefTetrahedron,order}) where order = 1
-Ferrite.nedgedofs(::Lagrange{2,RefTetrahedron,order}) where order = order - 1
+# Ferrite.nedgedofs(::Lagrange{2,RefTetrahedron,order}) where order = order - 1
Ferrite.nfacedofs(::Lagrange{2,RefTetrahedron,order}) where order = order - 1
+Ferrite.ncelldofs(::Lagrange{2,RefTetrahedron,3}) = 1
# permutation to switch to Ferrite ordering
permdof2D, permdof2Dinv = Dict{Int, Vector{Int}}(), Dict{Int, Vector{Int}}()Result: julia> test_H1Pk(degfem=3)
real.(λ) = [-7.083382139302583e-14, 2.467401100301959, 2.4674011003021055, 4.934802201505752, 9.869604408700198, 9.869604408705166, 12.337005525905012, 12.33700553941171, 19.739209047924067, 22.20661009725714]
nbdof = 8281
error = 3.8050011141876894e-8 |
|
Nice catch, thank you! I updated extend_interpolation.jl with the ordering required by Ferrite. The file |
|
.. I just read in @assert nfacedofs == 1 "Currently only supports interpolations with nfacedofs = 1"Is this limitation also for H^1 problems? |
| @@ -0,0 +1,83 @@ | |||
| function Ferrite.getnbasefunctions(::Lagrange{2,RefTetrahedron,order}) where order | |||
There was a problem hiding this comment.
If this file isn't used anymore, can you remove it?
src/interpolations.jl
Outdated
|
|
||
| #------------------------------------------------------------------------------- | ||
| # tentative to implement Lagrange element of order k ≥ 3 | ||
| include("extend_interpolations.jl") |
There was a problem hiding this comment.
I suggest you just add the content of this file directly after the current Lagrange{2,RefTetrahedron} in this file instead.
test/test_assembleP3.jl
Outdated
| return | ||
| end | ||
| #------------------------------------------------------------------------------- | ||
| function test_ip(ip) |
There was a problem hiding this comment.
I believe we have a test function like this in the test_interpolations.jl file, perhaps just test your new interpolations in that file?
|
I made some tests using another library to compare. I think the assembling part for order 4 and 5 may be incorrect. Perhaps due to |
|
What did you compare? I don't see anything wrong with this PR now, and dof distribution looks correct. |
Codecov ReportBase: 91.91% // Head: 90.75% // Decreases project coverage by
Additional details and impacted files@@ Coverage Diff @@
## master #482 +/- ##
==========================================
- Coverage 91.91% 90.75% -1.17%
==========================================
Files 22 22
Lines 3649 3699 +50
==========================================
+ Hits 3354 3357 +3
- Misses 295 342 +47
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report at Codecov. |
|
I used Getfem++ to assemble mass and stiffness matrices to solve the eigenvalue test problem we already discussed. |
|
Do you integrate it exactly ? The mass matrix with polynomial degree 4 would require Gauss integration of order 9, but we currently only have up to 8: At least when I run with order 8 for the quadrature I get: for example. |
|
Still the same mistake! |
|
I think order 8 is exact for degree 4. Order 10 would be required to be exact for order 5. |
IIRC, order |
|
You are right, I was thinking to the stiffness matrix |
|
I tried to follow the style of using PyCall, Printf
quadpy = pyimport("quadpy")
for order in 1:20
strelseif = order == 1 ? "if" : "elseif"
qpy = if order < 10
quadpy.t2.schemes["dunavant_0$(order)"]()
else
quadpy.t2.schemes["dunavant_$(order)"]()
end
nw = length(qpy.weights)
for k = 1:nw
if k == 1
println("$(strelseif) n == $(order)")
print("xw=[")
end
for l = 1:2
@printf("%1.16f ", qpy.points[l, k])
end
@printf("%1.16f / 2.0", qpy.weights[k])
k == nw ? println("];") : println("")
end
endwhich produces |
|
Good idea, perhaps that can be added in a separate PR. (I could not find the license for quadpy, but perhaps it doesn't matter when just using the output?) I will make some updates to this PR and then we can merge it I think. |
|
Hmm, not sure why GitHub closed the PR when I pushed some tests... I can't reopen it either. |
|
I think this is because you cannot push to remote master branches due to ownership conflicts. |
extend_interpolations.jlcontains the original version which does not exactly follow Ferrite ordering. The correspondence between reference_coordinates / values can be checked using thetest_ipfunction:test/test_assembleP3.jlalso contains an eigenvalue testtest_H1Pkwhich shows that theP3implementation is not workingextend_interpolations_Ferriteorder.jlis a tentative to reorder the basis to match Ferrite requirements. At this point it didn't solve the raised bytest_H1Pk