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

Interface elements #767

Closed
wants to merge 24 commits into from
Closed

Interface elements #767

wants to merge 24 commits into from

Conversation

DRollin
Copy link
Collaborator

@DRollin DRollin commented Jul 18, 2023

This PR implements cells, interpolations and cell values to be used for interface elements.
The main idea is that an interface element consists of two embedded elements.
This way the functionalities for embedded elements can be reused.
Since it does not affect the core functionalities and is more like a self-contained extension, I added the code in a separate file.

Some things still need some work in the future:

  • The interface should be consistent with what is being done for DG.
  • The way of referring to the two side needs to be discussed. Suggestions: :here&:there; :side1&:side2; true&false
  • A way of computing the normals of the interface will be needed for some applications and is not implemented yet.

src/interface_elements.jl Outdated Show resolved Hide resolved
src/interface_elements.jl Outdated Show resolved Hide resolved
src/interface_elements.jl Show resolved Hide resolved
src/interface_elements.jl Outdated Show resolved Hide resolved
src/interface_elements.jl Outdated Show resolved Hide resolved
Co-authored-by: Fredrik Ekre <[email protected]>
@codecov-commenter
Copy link

codecov-commenter commented Jul 20, 2023

Codecov Report

Attention: Patch coverage is 99.00498% with 2 lines in your changes are missing coverage. Please review.

Project coverage is 93.53%. Comparing base (8e92196) to head (71914da).
Report is 1 commits behind head on master.

Files Patch % Lines
src/interface_elements.jl 99.00% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master     #767      +/-   ##
==========================================
+ Coverage   93.30%   93.53%   +0.22%     
==========================================
  Files          36       37       +1     
  Lines        5257     5459     +202     
==========================================
+ Hits         4905     5106     +201     
- Misses        352      353       +1     

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

DRollin and others added 3 commits July 20, 2023 07:56
Comment on lines 85 to 95
function get_node_ids(c::InterfaceCell)
return collect( getproperty(c, side).nodes[baseindex] for (side, baseindex) in get_sides_and_base_indices(c) )
end

function Base.getproperty(c::InterfaceCell, s::Symbol)
if s == :nodes
return (c.here.nodes..., c.there.nodes...) # Order used to reinit InterfaceCellValues
else
return getfield(c,s)
end
end
Copy link
Member

Choose a reason for hiding this comment

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

Can you comment on this. These seem to return in different orders, right? I guess it is intentional, but a comment explaining why would be nice.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I wanted to avoid that the correct node indices have to be figured out at every reinit!-call and was a bit lazy with the solution.
To be consistent with the usual ordering, I changed this:
Now the required indices for sorting the coordinates are determined when creating the InterfaceCellValues and stored in vectors.
I am not sure if there is a better way of doing it, but this way, this is at least only done once.

Copy link
Collaborator

@lijas lijas left a comment

Choose a reason for hiding this comment

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

Nice work!

I guess functions like function_jump_value and function_jump_gradient could/should also be added?

src/interface_elements.jl Outdated Show resolved Hide resolved
src/interface_elements.jl Outdated Show resolved Hide resolved
@kimauth
Copy link
Member

kimauth commented Mar 27, 2024

I can't distribute dofs on 2D interface elements:

dim = 2
grid = generate_grid(Quadrilateral, (3,1))

# Ferrite default grid:
# 5 ____ 6 ____ 7 ____ 8
# |      |      |      |
# |      |      |      |
# 1 ____ 2 ____ 3 ____ 4

# make middle cell zero-thickness:
grid.nodes[2:3] = [Node(mean(n->n.x, grid.nodes[2:3])) for _ in 1:2]
grid.nodes[6:7] = [Node(mean(n->n.x, grid.nodes[6:7])) for _ in 1:2]

# 1D-line elements for construction of InterfaceCell
line1 = Line((2,6))
line2 = Line((3,7))

cohesivecell = InterfaceCell(line1, line2)

cohesive_grid = Grid([grid.cells[1], grid.cells[3], cohesivecell], grid.nodes)

ip_quads = Lagrange{RefQuadrilateral, 1}()#^dim
ip_base = Lagrange{RefLine, 1}()
ip_cohesive = InterfaceCellInterpolation(ip_base)#^dim

dh = DofHandler(cohesive_grid)
sdh_quads = SubDofHandler(dh, Set(1:2))
add!(sdh_quads, :u, ip_quads)
sdh_coh = SubDofHandler(dh, Set(3))
add!(sdh_coh, :u, ip_cohesive)
close!(dh)
julia> close!(dh)                                                                                                                                                                                                                                       [3/1802]
ERROR: MethodError: no method matching facedof_interior_indices(::Lagrange{RefLine, 1, Nothing})                                                                                                                                                                
                                                                                                                                                                                                                                                                
Closest candidates are:                                                                                                                                                                                                                                         
  facedof_interior_indices(::InterfaceCellInterpolation)
   @ Ferrite ~/.julia/packages/Ferrite/9NWtq/src/interface_elements.jl:215
  facedof_interior_indices(::CrouzeixRaviart)
   @ Ferrite ~/.julia/packages/Ferrite/9NWtq/src/interpolations.jl:1296
  facedof_interior_indices(::Lagrange{RefQuadrilateral, 2})
   @ Ferrite ~/.julia/packages/Ferrite/9NWtq/src/interpolations.jl:530

"""
function get_interface_index(ip::InterfaceCellInterpolation, side::Symbol, i::Integer)
nvhere, nvthere = get_n_dofs_on_side(vertexdof_indices, ip, :here), get_n_dofs_on_side(vertexdof_indices, ip, :there)
nfhere, nfthere = get_n_dofs_on_side(facedof_interior_indices, ip, :here), get_n_dofs_on_side(facedof_interior_indices, ip, :there)
Copy link
Member

Choose a reason for hiding this comment

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

I think this needs some kind of branch for dim=2

Copy link
Member

Choose a reason for hiding this comment

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

I think we should consider defining empty defaults such as

Ferrite.facedof_interior_indices(::Lagrange{RefLine,1}) = ((), ())

instead of introducing branches. This would just mean that a line does not have interior facedofs. Imo that's easier to understand than branches.

(Introducing this fixes the problem for 2d dof distribution reported above).

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I agree 👍
Would
Ferrite.facedof_interior_indices(::Lagrange{RefLine}) = ((), ())
be more general?


vertexdof_indices(ip::InterfaceCellInterpolation) = get_interface_dof_indices(vertexdof_indices, ip)

function facedof_indices(ip::InterfaceCellInterpolation)
Copy link
Member

Choose a reason for hiding this comment

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

Can you clarify what face means in this context?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I think of the two sides of the interface as faces.
In 2D that means the interface is like a Quadrilateral, but without two faces and no interior due to the zero thickness of the element.

src/exports.jl Outdated Show resolved Hide resolved
InterfaceCell,
InterfaceCellInterpolation,
InterfaceCellValues,
get_side_and_baseindex,
Copy link
Member

Choose a reason for hiding this comment

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

Is this user facing?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Except for get_side_and_baseindex definitely.
That function is something I planned to use in the future, where I need to evaluate values on both sides separately.
This is probably a quite specific case, but I thought that I might not be the only one using it like this at some point.

src/interface_elements.jl Outdated Show resolved Hide resolved
cell = InterfaceCell(here, there)

@test Ferrite.nvertices(cell) == Ferrite.nvertices(here) + Ferrite.nvertices(there)
@test Ferrite.nfaces(cell) == 2
Copy link
Member

Choose a reason for hiding this comment

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

Where does this assumption come from?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

The idea is that the two sides of an interface are considered the faces of the cell.
To me it makes sense to think of it this way, but theses definitions can definitely be discussed.

@DRollin
Copy link
Collaborator Author

DRollin commented Apr 19, 2024

Thank you all for your feedback!
After some discussions, we decided to create a separate package for the interface elements: FerriteInterfaceElements.jl
Since it is an add-on and not a core functionality, we thought this a better way to do it.
Hence, I will close this PR.

@DRollin DRollin closed this Apr 19, 2024
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.

7 participants