Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions src/Dofs/DofHandler.jl
Original file line number Diff line number Diff line change
Expand Up @@ -644,7 +644,7 @@ described therein.
# typestable way.
n_copies = step(dofs)
@assert n_copies > 0
if adjust_during_distribution && !orientation.regular
if adjust_during_distribution && orientation.flipped
# Reverse the dofs for the path
dofs = reverse(dofs)
end
Expand All @@ -661,12 +661,12 @@ end

Returns the unique representation of an edge and its orientation.
Here the unique representation is the sorted node index tuple. The
orientation is `true` if the edge is not flipped, where it is `false`
orientation is `false` if the edge is not flipped, where it is `true`
if the edge is flipped.
"""
function sortedge(edge::Tuple{Int, Int})
a, b = edge
return a < b ? (edge, PathOrientationInfo(true)) : ((b, a), PathOrientationInfo(false))
return a < b ? (edge, PathOrientationInfo(false)) : ((b, a), PathOrientationInfo(true))
end

"""
Expand Down Expand Up @@ -739,7 +739,7 @@ function sortface(face::Tuple{Int, Int, Int})
b, c = minmax(b, c)
a, c = minmax(a, c)
a, b = minmax(a, b)
return (a, b, c), SurfaceOrientationInfo() # TODO fill struct
return (a, b, c), SurfaceOrientationInfo(face)
end


Expand All @@ -760,7 +760,7 @@ function sortface(face::Tuple{Int, Int, Int, Int})
b, c = minmax(b, c)
a, c = minmax(a, c)
a, b = minmax(a, b)
return (a, b, c), SurfaceOrientationInfo() # TODO fill struct
return (a, b, c), SurfaceOrientationInfo(face)
end


Expand Down
4 changes: 2 additions & 2 deletions src/FEValues/InterfaceValues.jl
Original file line number Diff line number Diff line change
Expand Up @@ -419,8 +419,8 @@ Return the relative orientation info for facet B with regards to facet A.
Relative orientation is computed using a [`OrientationInfo`](@ref) for each side of the interface.
"""
function InterfaceOrientationInfo(cell_a::AbstractCell{RefShapeA}, cell_b::AbstractCell{RefShapeB}, facet_a::Int, facet_b::Int) where {RefShapeA <: AbstractRefShape, RefShapeB <: AbstractRefShape}
OI_a = OrientationInfo(facets(cell_a)[facet_a])
OI_b = OrientationInfo(facets(cell_b)[facet_b])
OI_a = orientation_info(facets(cell_a)[facet_a])
OI_b = orientation_info(facets(cell_b)[facet_b])
flipped = OI_a.flipped != OI_b.flipped
shift_index = OI_b.shift_index - OI_a.shift_index
return InterfaceOrientationInfo{RefShapeA, RefShapeB}(flipped, shift_index, OI_b.shift_index, facet_a, facet_b)
Expand Down
126 changes: 51 additions & 75 deletions src/Grid/grid.jl
Original file line number Diff line number Diff line change
Expand Up @@ -664,8 +664,28 @@ end
#################################
# @TODO merge this code with into the logic in `ConstraintHandler`.

function get_edge_direction(edgenodes::NTuple{2, Int})
positive = edgenodes[2] > edgenodes[1]
return ifelse(positive, 1, -1)
end

function get_face_direction(facenodes::NTuple{N, Int}) where {N}
N > 2 || throw(ArgumentError("A face must have at least 3 nodes"))
min_idx = argmin(facenodes)
if min_idx == 1
positive = facenodes[2] < facenodes[end]
elseif min_idx == length(facenodes)
positive = facenodes[1] < facenodes[end - 1]
else
positive = facenodes[min_idx + 1] < facenodes[min_idx - 1]
end
return ifelse(positive, 1, -1)
end

"""
PathOrientationInfo
PathOrientationInfo(flipped:Bool)
PathOrientationInfo(edgenodes::NTuple{2, Int})

Orientation information for 1D entities.

Expand All @@ -674,117 +694,73 @@ associated to the vertices. To give an example, the oriented path
```
1 ---> 2
```
is called *regular*, indicated by `regular=true`, while the oriented path
is called *regular*, indicated by `flipped=false`, while the oriented path
```
2 ---> 1
```
is called *inverted*, indicated by `regular=false`.
is called *inverted*, indicated by `flipped=true`.
"""
struct PathOrientationInfo
regular::Bool # Indicator whether the orientation is regular or inverted.
flipped::Bool # Indicator whether the orientation is regular or inverted.
Copy link
Member

Choose a reason for hiding this comment

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

IMO it is less straight forward to store if the edge is flipped, than if it is regular. Personally, I would prefer positive, or perhaps we even can save simply as

direction::Int8 # \pm 1

and return Int8 from get_face_direction and get_edge_direction. Probably, we won't see any significant performance difference by using Int8 vs Int, don't recall the profiles of the dof distribution now if this logic takes some significant portion of the time...

Copy link
Member

Choose a reason for hiding this comment

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

Reagardng the dof distribution time, it is mostly spend in the hash table lookups, as we do not materialize the mesh entities explicitly. But I think just a boolean is clearer here. Also, I would expect that flipped = true would mean 2--->1.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

A "flip" is indeed the terminology used in the Scroggs paper to indicate that the nodes are swapped compared to the "normal" orientation. So,

  • 1 --> 2 is flipped = false
  • 2 --> 1 is flipped = true

This is the same way OrientationInfo already worked and in line with the paper. The previous definition of regular::Bool was opposite, of course. Can you check whether the tests make sense following this definition?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I see now that the docstring was not updated yet.
Refer to the latest commit for the corrected definition.

shift_index::Int

function PathOrientationInfo(flipped::Bool)
return new(flipped, 0)
end
end

function PathOrientationInfo(edgenodes::NTuple{2, Int})
return PathOrientationInfo(get_edge_direction(edgenodes) < 0)
end

"""
SurfaceOrientationInfo
SurfaceOrientationInfo(flipped:Bool, shift_index::Int)
SurfaceOrientationInfo(facenodes::NTuple{N, Int})

Orientation information for 2D entities. Such an entity can be
possibly flipped (i.e. the defining vertex order is reverse to the
spanning vertex order) and the vertices can be rotated against each other.
Take for example the faces
```
1---2 2---3
4---3 3---2
| A | | B |
4---3 1---4
1---2 4---1
```
which are rotated against each other by 90° (shift index is 1) or the faces
```
1---2 2---1
4---3 2---3
| A | | B |
4---3 3---4
1---2 1---4
```
which are flipped against each other. Any combination of these can happen.
The combination to map this local face to the defining face is encoded with
this data structure via ``rotate \\circ flip`` where the rotation is indiced by
the shift index.
!!!NOTE TODO implement me.
"""
struct SurfaceOrientationInfo
#flipped::Bool
#shift_index::Int
end


@doc raw"""
InterfaceOrientationInfo

Orientation information for 1D and 2D entities.
The orientation is defined by the indices of the grid nodes
associated to the vertices. To give an example, the oriented path
```
1 ---> 2
```
is called *regular*, indicated by `flipped=false`, while the oriented path
```
2 ---> 1
```
is called *inverted*, indicated by `flipped=true`.

2D entities can be flipped (i.e. the defining vertex order is reverse to the
spanning vertex order) and the vertices can be rotated against each other.

The reference entity is a one with it's first node is the lowest index vertex
and its vertices span counter-clock-wise.
Take for example the faces
```
1 2
| \ | \
| \ | \
| A \ | B \
| \ | \
2-----3 3-----1
```
which are rotated against each other by 240° after tranfroming to an
equilateral triangle (shift index is 2). Or the faces
```
3 2
| \ | \
| \ | \
| A \ | B \
| \ | \
2-----1 3-----1
```
which are flipped against each other.
"""
struct OrientationInfo
flipped::Bool
shift_index::Int
end

function OrientationInfo(edgenodes::NTuple{2, Int})
return OrientationInfo(get_edge_direction(edgenodes) < 0, 0)
end

function OrientationInfo(facenodes::NTuple{N, Int}) where {N}
function SurfaceOrientationInfo(facenodes::NTuple{N, Int}) where {N}
min_idx = argmin(facenodes)
shift_index = min_idx - 1
flipped = get_face_direction(facenodes) < 0
return OrientationInfo(flipped, shift_index)
return SurfaceOrientationInfo(flipped, shift_index)
end

function get_edge_direction(edgenodes::NTuple{2, Int})
positive = edgenodes[2] > edgenodes[1]
return ifelse(positive, 1, -1)

"""
orientation_info(edgenodes::NTuple{2, Int})
orientation_info(facenodes::NTuple{N, Int})

Helper function that returns a `PathOrientationInfo` or `SurfaceOrientationInfo` as appropriate for the size of the node tuple.
"""
function orientation_info(edgenodes::NTuple{2, Int})
return PathOrientationInfo(edgenodes)
end

function get_face_direction(facenodes::NTuple{N, Int}) where {N}
N > 2 || throw(ArgumentError("A face must have at least 3 nodes"))
min_idx = argmin(facenodes)
if min_idx == 1
positive = facenodes[2] < facenodes[end]
elseif min_idx == length(facenodes)
positive = facenodes[1] < facenodes[end - 1]
else
positive = facenodes[min_idx + 1] < facenodes[min_idx - 1]
end
return ifelse(positive, 1, -1)
function orientation_info(facenodes::NTuple{N, Int}) where {N}
return SurfaceOrientationInfo(facenodes)
end
1 change: 1 addition & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ include("test_utils.jl")
# Unit tests
include("test_collectionsofviews.jl")
include("test_refshapes.jl")
include("test_orientations.jl")
include("test_interpolations.jl")
include("test_cellvalues.jl")
include("test_facevalues.jl")
Expand Down
75 changes: 75 additions & 0 deletions test/test_orientations.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
@testset "Orientations" begin
@testset "PathOrientationInfo" begin
# Path [a, b], flipped
test_paths = [
((1, 2), false), # 1 -> 2, regular
((2, 1), true), # 2 -> 1, inverted
((5, 1), true), # 5 -> 1, inverted
((3, 5), false), # 3 -> 5, regular
]

for (path, flipped) in test_paths
orientation = Ferrite.PathOrientationInfo(path)

@test orientation.flipped == flipped
end
end

@testset "SurfaceOrientationInfo" begin
# Surface [facenodes], flipped, shift_index
test_surfaces = [
# 2
# | \
# 3 - 1, regular, rotation = 0 deg
((1, 2, 3), false, 0),
# 1
# | \
# 2 - 3, regular, rotation = 120 deg
((3, 1, 2), false, 1),
# 3
# | \
# 1 - 2, regular, rotation = 240 deg
((2, 3, 1), false, 2),
# 3
# | \
# 2 - 1, inverted, rotation = 0 deg
((1, 3, 2), true, 0),
# 1
# | \
# 3 - 2, inverted, rotation = 120 deg
((2, 1, 3), true, 1),
# 2
# | \
# 1 - 3, inverted, rotation = 240 deg
((3, 2, 1), true, 2),
# 4 - 3
# | |
# 1 - 2, regular, rotation = 0 deg
((1, 2, 3, 4), false, 0),
# 2 - 3
# | |
# 1 - 4, inverted, rotation = 0 deg
((1, 4, 3, 2), true, 0),
# 3 - 2
# | |
# 4 - 1, regular, rotation = 90 deg
((4, 1, 2, 3), false, 1),
# 1 - 2
# | |
# 4 - 3, inverted, rotation = 270 deg
((4, 3, 2, 1), true, 3),
]


for (surface, flipped, shift_index) in test_surfaces
orientation = Ferrite.SurfaceOrientationInfo(surface)

@test orientation.flipped == flipped
@test orientation.shift_index == shift_index
end
end

@assert typeof(Ferrite.orientation_info((1, 2))) <: Ferrite.PathOrientationInfo
@assert typeof(Ferrite.orientation_info((1, 2, 3))) <: Ferrite.SurfaceOrientationInfo
@assert typeof(Ferrite.orientation_info((1, 2, 3, 4))) <: Ferrite.SurfaceOrientationInfo
end
Loading