Skip to content
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

Wedge Support #581

Merged
merged 46 commits into from
Mar 17, 2023
Merged

Wedge Support #581

merged 46 commits into from
Mar 17, 2023

Conversation

termi-official
Copy link
Member

@termi-official termi-official commented Jan 10, 2023

This PR contains a fully functional code to support linear and quadratic wedges (see #580) for both dof handlers. Addresses #272 as a side-effect.

@termi-official
Copy link
Member Author

Btw the fail in the CI is unrelated to this PR, but I found it when improving the tests. The way the quintic triangle is evaluated is not really stable (numerically).

@KristofferC
Copy link
Collaborator

Cool, this is an important case for improved generality.

@termi-official
Copy link
Member Author

termi-official commented Jan 12, 2023

CI is broken due to the deprecation. Refinement leads to decrease in the L_infty norm on Poisson problems with manufactured solutions, indicating convergence.

Edit: Test coverage for the new elements is incomplete right now and will be completed when the review progresses.

Edit 2: Linear wedges are exportable to VTK. To visually inspect the quadratic wedges please use FerriteViz#do/crincle-clip (Ferrite-FEM/FerriteViz.jl#56)

together with this code to make wedges functional

function FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Ferrite.AbstractCell{3,N,5}) where {N}
    faces = Ferrite.faces(cell)
    (coord_offset, triangle_offset) = FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, Ferrite.Cell{3,3,1}(faces[1]))
    (coord_offset, triangle_offset) = FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, Ferrite.Cell{3,4,1}(faces[2]))
    (coord_offset, triangle_offset) = FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, Ferrite.Cell{3,4,1}(faces[3]))
    (coord_offset, triangle_offset) = FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, Ferrite.Cell{3,4,1}(faces[4]))
    (coord_offset, triangle_offset) = FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, Ferrite.Cell{3,3,1}(faces[5]))
    (coord_offset, triangle_offset)
end

FerriteViz.ntriangles(cell::Ferrite.AbstractCell{3,N,5}) where {N} = 1+4+4+4+1

FerriteViz.getfaceip(ip::Ferrite.Lagrange{3, RefPrism, order}, local_face_idx::Int) where {order} = local_face_idx in 2:4 ? Ferrite.Lagrange{2, RefCube, order}() : Ferrite.Lagrange{2, RefTetrahedron, order}()

FerriteViz.linear_face_cell(cell::Ferrite.Cell{3,N,5}, local_face_idx::Int) where N = local_face_idx in 2:4 ? Ferrite.Quadrilateral3D(Ferrite.faces(cell)[local_face_idx]) : Ferrite.Cell{3,3,1}(Ferrite.faces(cell)[local_face_idx])

Copy link
Member

@KnutAM KnutAM left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cool!
I looked through, out of interest, and spotted some minor things.

src/FEValues/face_values.jl Outdated Show resolved Hide resolved
src/Dofs/DofHandler.jl Outdated Show resolved Hide resolved
src/interpolations.jl Outdated Show resolved Hide resolved
@lijas lijas mentioned this pull request Feb 28, 2023
@termi-official
Copy link
Member Author

The faulty assumption

edges(ip::Interpolation{2}) = faces(ip)
edgedof_indices(ip::Interpolation{2}) = facedof_indices(ip)
edgedof_interior_indices(ip::Interpolation{2}) = facedof_interior_indices(ip)

has to be fixed in a subsequent PR. I could not find a simple way to propagate the shell dof elimination through the constraint handler without breaking dof elimination for mixed grids.

@codecov-commenter
Copy link

codecov-commenter commented Mar 3, 2023

Codecov Report

Patch coverage: 81.32% and project coverage change: -0.02 ⚠️

Comparison is base (200b47a) 92.66% compared to head (70a9c8a) 92.65%.

📣 This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #581      +/-   ##
==========================================
- Coverage   92.66%   92.65%   -0.02%     
==========================================
  Files          28       29       +1     
  Lines        4270     4441     +171     
==========================================
+ Hits         3957     4115     +158     
- Misses        313      326      +13     
Impacted Files Coverage Δ
src/Export/VTK.jl 89.83% <0.00%> (-1.55%) ⬇️
src/Quadrature/gaussquad_prism_table.jl 0.00% <0.00%> (ø)
src/Quadrature/quadrature.jl 77.61% <0.00%> (-16.94%) ⬇️
src/Dofs/MixedDofHandler.jl 84.82% <83.78%> (-0.83%) ⬇️
src/Dofs/DofHandler.jl 89.80% <90.47%> (-0.66%) ⬇️
src/interpolations.jl 97.90% <93.75%> (+7.43%) ⬆️
src/Dofs/ConstraintHandler.jl 96.05% <100.00%> (+1.11%) ⬆️
src/FEValues/face_values.jl 100.00% <100.00%> (ø)
src/Ferrite.jl 100.00% <100.00%> (ø)
src/Grid/grid.jl 90.88% <100.00%> (ø)
... and 2 more

... and 1 file with indirect coverage changes

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 in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

@lijas
Copy link
Collaborator

lijas commented Mar 3, 2023

I think #477 handles shell interpolations better. Maybe merging that PR will makes life easier in this PR?

@termi-official
Copy link
Member Author

It would definitely make things easier, but the underlying issue that the semantics of edge and face is different in different sections of the codebase stays.
CI works again, so this PR should be ready.

Copy link
Member Author

@termi-official termi-official left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should summarize the remaining issues and clarify some stuff @fredrikekre

docs/src/devdocs/interpolations.md Outdated Show resolved Hide resolved
src/Dofs/DofHandler.jl Show resolved Hide resolved
src/Dofs/DofHandler.jl Show resolved Hide resolved
src/Dofs/MixedDofHandler.jl Show resolved Hide resolved
src/FEValues/face_values.jl Show resolved Hide resolved
@@ -154,7 +179,7 @@ indices of [`reference_coordinates(::Interpolation)`](@ref).
value(ip::Interpolation, i::Int, ξ::Vec)

"""
reference_coordinates(::Interpolation)
Ferrite.reference_coordinates(::Interpolation)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When should we include the Ferrite. prefix in the doc string?

TODO for me: Put answer into devdocs.

src/interpolations.jl Outdated Show resolved Hide resolved
src/interpolations.jl Outdated Show resolved Hide resolved
test/test_interpolations.jl Show resolved Hide resolved
test/test_interpolations.jl Show resolved Hide resolved
@fredrikekre fredrikekre merged commit 1c01761 into master Mar 17, 2023
@fredrikekre fredrikekre deleted the do/wedge-elements branch March 17, 2023 16:11
@termi-official
Copy link
Member Author

Thanks for working through this !

@fredrikekre fredrikekre restored the do/wedge-elements branch March 17, 2023 22:01
@fredrikekre fredrikekre deleted the do/wedge-elements branch March 17, 2023 22:06
fredrikekre added a commit that referenced this pull request Apr 18, 2023
This patch reworks the `AbstractCell` interface (refer to #677 for
motivation and discussion for alternatives to this patch). Specifically,
this patch

 - introduces the reference dimension as a type parameter on `AbstractRefShape`,
 - introduces `RefInterval`, `RefTriangle`, `RefQuadrilateral`, and
   `RefHexahedron` as new subtypes of `AbstractRefShape{refdim}` (note
   that interpolations are not updated to this system yet since it can
   be done in a separate patch),
 - deprecates the old `Cell{refdim, nnodes, nfaces}` struct,
 - introduces new structs for all supported cell types, e.g.
   `struct Triangle <: AbstractCell{RefTriangle}`,
 - defines the unique vertex/edge/face identifiers for DoF-distribution
   in terms of the reference shape.

The new parametrization makes it possible to dispatch on the reference
shape instead of union of concrete implementations, i.e.
`::AbstractCell{RefTriangle}` instead of `::Union{Triangle,
QuadraticTriangle}`, and similarly for the reference dimension.

One drawback is that it is no longer possible to use "anonymous
celltypes" after this patch, i.e. using, e.g., `Cell{3, 20, 6}` and
defining methods for that instead of explicitly naming it
`SerendipityQuadraticHexahedron` or similar. Instead, you now have to
define a struct and some methods. Arguably this is a good thing --
nobody understands the meaning of `Cell{3, 20, 6}`.

For normal usage this patch is not breaking (see e.g. no required
changes for examples), unless one dispatches on the `Cell` struct
directly.

Specific things that remain to be decided:

 - `struct Triangle <: AbstractCell{2, RefTriangle}` or `struct Triangle
   <: AbstractRefCell{2, RefTriangle} <: AbstractCell`? I don't know if
   it is too restrictive to "force" `{refdim, refshape}` on all subtypes
   of `AbstractCell`. On the other hand, if one implements a
   non-refshape based element one can define a dummy reference shape.
   For a real-world example, looking at e.g.
   https://github.com/kimauth/FerriteCohesiveZones.jl/blob/main/src/cohesive_cell.jl
   I think that would simply be implemented as
   ```julia
   struct CohesiveQuadrilateral <: AbstractCell{2, RefQuadrilateral}
       nodes::NTuple{4, Int}
   end
   ```
   which I don't see a direct problem with, other than the fact that the
   cell will now have 4 edges (according to the fallback definitions).
   To solve this, one would have to implement a method for
   `faces(::CohesiveQuadrilateral)` that only return the upper and lower
   edge. I don't think this is too much to ask for custom elements like
   this. Alternatively, with #581 merged it is possible to distribute
   different number of dofs for different edges, so this can also be
   solved on the (default) interpolation level by distributing 0 dofs on
   the left and right edge.

 - `AbstractCell{refdim, refshape}` or just `AbstractCell{refshape}`?
   The dimension is implicit from the shape now, and you can dispatch on
   `AbstractCell{<:AbstractRefShape{refdim}} where refdim`. Personally I
   am in favor of including the dimension, even if redundant. It does
   not cost anything, and sometimes you want to dispatch on the
   dimension, sometimes on the refshape.

Fixes #677.
fredrikekre added a commit that referenced this pull request Apr 24, 2023
After #581 `vertexdof_indices` should return a tuple of tuples, but
since `length(::Int) == length(::Tuple{Int}) == 1` this old code worked
anyway. This patch removes the method to rely on the generic `Lagrange`
method instead.
fredrikekre added a commit that referenced this pull request Apr 24, 2023
After #581 `vertexdof_indices` should return a tuple of tuples, but
since `length(::Int) == length(::Tuple{Int}) == 1` this old code worked
anyway. This patch removes the method to rely on the generic `Lagrange`
method instead.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants