check_isoparametric_boundaries, mass_qr, and cellcenter for Pyramids and Wedges + L2Projector tests for Hexahedron, Tetrahedron, Pyramid, and Wedge.#1215
Open
AbdAlazezAhmed wants to merge 10 commits intoFerrite-FEM:masterfrom
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1215 +/- ##
==========================================
- Coverage 94.25% 93.07% -1.18%
==========================================
Files 40 40
Lines 6750 6859 +109
==========================================
+ Hits 6362 6384 +22
- Misses 388 475 +87 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
| [deps] | ||
| EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" | ||
| ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" | ||
| Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" |
lijas
reviewed
Aug 12, 2025
Collaborator
lijas
left a comment
lijas
left a comment
There was a problem hiding this comment.
Is it easy to implement the linear interpolation ourself and avoid importing the Interpolations package?
| function test_projection(order, elementtype) | ||
| refshape = getrefshape(elementtype) | ||
| dim = Ferrite.getrefdim(refshape) | ||
| grid = single_element_grid(elementtype) |
Collaborator
There was a problem hiding this comment.
Why did you change this to a single element grid? Previously we had triangles/tetrahedrons filled a square/cube.
Collaborator
Author
There was a problem hiding this comment.
My bad, I didn't notice it.
| test_projection(2, RefQuadrilateral) | ||
| test_projection(2, RefTriangle) | ||
| for ref_shape in ( | ||
| Quadrilateral, |
| # A linear approximation can not recover a quadratic solution, | ||
| # so projected values will be different from the analytical ones | ||
| ae = [f_approx(i) for i in 1:4] | ||
| qp_1D_coord = Float64[] |
Collaborator
There was a problem hiding this comment.
Add a small comment what you are trying to do here
Collaborator
Author
There was a problem hiding this comment.
I was trying to calculate the first order approximation values using linear extrapolation from the quadrature points instead of using precomputed hardcoded values. But since we use more than one element this won't work so I'll compute the values for other cell types and add them instead.
…mids and prisms need to be increased so that the matrix is matrixing
AbdAlazezAhmed
changed the title
check_isoparametric_boundaries and cellcenter for Pyramids and Wedges + better tests for L2Projector.check_isoparametric_boundaries, mass_qr, and cellcenter for Pyramids and Wedges + L2Projector tests for Hexahedron, Tetrahedron, Pyramid, and Wedge.
Member
The new functions are never called in the tests. |
Collaborator
Author
|
So it seems like there's something off with the quadrature rules for prisms and pyramids. Here's a copied test to check what quadrature orders work: Some Copied Codeusing Ferrite
function test_pe_scalar_field(celltype::Type{<:Ferrite.AbstractCell}, ip_order, qr_order)
ref_shape = getrefshape(celltype)
dim = Ferrite.getrefdim(ref_shape)
mesh = generate_grid(celltype, ntuple(_ -> 3, dim))
f(x) = sum(x)
ip_f = Lagrange{ref_shape, ip_order}() # function interpolation
ip_g = geometric_interpolation(celltype) # geometry interpolation
# points where we want to retrieve field values
points = Vec{dim, Float64}[]
# compute values in quadrature points
qr = Ferrite.QuadratureRule{ref_shape}(qr_order) # exactly integrate field
cv = CellValues(qr, ip_f, ip_g)
qp_vals = [Vector{Float64}(undef, getnquadpoints(cv)) for _ in 1:getncells(mesh)]
for cellid in eachindex(mesh.cells)
xe = getcoordinates(mesh, cellid)
reinit!(cv, xe)
for qp in 1:getnquadpoints(cv)
x = spatial_coordinate(cv, qp, xe)
qp_vals[cellid][qp] = f(x)
push!(points, x)
end
end
# do a L2Projection for getting values in dofs
projector = L2Projector(ip_f, mesh)
projector_vals = project(projector, qp_vals, qr)
return maximum(projector_vals) ≈ dim
end
for ip_order in 1:3
for qr_order in 1:7
function Ferrite._mass_qr(::Lagrange{RefPyramid, order}) where {order}
return QuadratureRule{RefPyramid}(qr_order)
end
function Ferrite._mass_qr(::Lagrange{RefPrism, order}) where {order}
return QuadratureRule{RefPrism}(qr_order)
end
for cell_type in (Wedge, Pyramid)
try
worked = test_pe_scalar_field(cell_type, ip_order, qr_order)
println("cell: $cell_type, ip order: $ip_order, qr order: $qr_order ::: correct values? $worked")
catch
println("cell: $cell_type, ip order: $ip_order, qr order: $qr_order ::: Crashed")
end
println("--------------")
end
end
end
Results: All results |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
PointEvalHandler needed these functions to work. Also, L2 Projector is only tested with triangles and quads. I added a bit more general test using
Interpolations.jlinstead of manually calculating the test values for first order interpolations.