Open fredrikekre opened 3 months ago
I could probably have a go at this myself with a bit of help understanding the Grid
struct. I may be the wrong person for it since I've never actually used this package, but maybe it's useful. A realistic example could be a curve-bounded domain, e.g. https://juliageometry.github.io/DelaunayTriangulation.jl/stable/tutorials/curve_bounded/#A-Complicated-Multiply-Connected-Disjoint-Domain (not this domain exactly since I don't imagine it's very practical, but just to show what could be done).
For example of what I'm having difficulty understanding, with
tri = triangulate_rectangle(0, 1, 0, 1, 4, 4, single_boundary = false) # [0, 1] × [0, 1], 20 points in each direction
triplot(tri)
each_solid_triangle(tri) |> collect
18-element Vector{Tuple{Int64, Int64, Int64}}:
(9, 10, 13)
(1, 2, 5)
(7, 4, 8)
(6, 3, 7)
(7, 8, 11)
(11, 8, 12)
(9, 6, 10)
(2, 3, 6)
(10, 7, 11)
(11, 12, 15)
(3, 4, 7)
(5, 2, 6)
(6, 7, 10)
(13, 10, 14)
(5, 6, 9)
(14, 11, 15)
(10, 11, 14)
(15, 12, 16)
The cells can be easily obtained by iterating over each_solid_triangle(tri)
, and the nodes from each_solid_vertex
and get_point
. (If you are unfamiliar, the corresponding _ghost
iterators are used for boundary information, as used in e.g. FiniteVolumeMethod.jl.) I don't really understand what (:cellsets, :nodesets, :facetsets, :vertexsets)
are defined as exactly, nor their ordering (which I imagine is important given the use of OrderedSets
). For example,
grid = generate_grid(Triangle, (4, 4))
julia> grid.cellsets
Dict{String, OrderedCollections.OrderedSet{Int64}}()
julia> grid.nodesets
Dict{String, OrderedCollections.OrderedSet{Int64}}()
julia> grid.facetsets
Dict{String, OrderedCollections.OrderedSet{FacetIndex}} with 4 entries:
"left" => OrderedSet{FacetIndex}([FacetIndex((1, 3)), FacetIndex((9, 3)), FacetIndex((17, 3)), FacetIndex((25, 3))])
"bottom" => OrderedSet{FacetIndex}([FacetIndex((1, 1)), FacetIndex((3, 1)), FacetIndex((5, 1)), FacetIndex((7, 1))])
"right" => OrderedSet{FacetIndex}([FacetIndex((8, 1)), FacetIndex((16, 1)), FacetIndex((24, 1)), FacetIndex((32, 1))])
"top" => OrderedSet{FacetIndex}([FacetIndex((26, 2)), FacetIndex((28, 2)), FacetIndex((30, 2)), FacetIndex((32, 2))])
julia> grid.vertexsets
Dict{String, OrderedCollections.OrderedSet{VertexIndex}}()
I'm not sure what I'm actually looking at here. I did see https://ferrite-fem.github.io/Ferrite.jl/dev/topics/grid/, but I don't really follow it unfortunately.
I'm also assuming that the aim would be to put all the information into Grid
at least, which is probably simplest just for a one-off example. Another way would be to define a wrapper
struct DelaunayGrid{T<:Triangulation} <: AbstractGrid
triangulation::T
end
and use that as the interface. Could be for a useful package extension in the future maybe.
Cool!
I don't really understand what (:cellsets, :nodesets, :facetsets, :vertexsets) are defined as exactly
These sets define subsets of the grid, to allow performing operations on only the subset.
The two most common to be filled are cellsets
and facetsets
. cellsets
simply gives a list of cell numbers that can be used to perform specific operations on part of the domain (for example having different physical properties). For making this example, I don't think that will be necessary.
The facetsets
are typically used to define boundary conditions for the simulation. The FacetIndex
is basically a tuple of (cell_number, local_facet_number), where in 2d the local_facet_number
implies which of the edges of the cell the FacetIndex
refers to.
However, you do not have to populate this either, since that can be done by calling add_facetset!
or add_boundary_facetset!
after having created the grid. With your boundary iterator it might be possible to make more efficient, but that could be added later!
The ordering of any of these should have no effect on the simulation results, but may affect performance during assembly and the resulting finite element matrix structure (for cellset).
The FacetIndex is basically a tuple of (cell_number, local_facet_number), where in 2d the local_facet_number implies which of the edges of the cell the FacetIndex refers to.
I see. In grid.facetsets
in my example, I see that the left
boundary has a FacetIndex(1, 3)
, so
julia> cn, lcn = FacetIndex(1, 3)
FacetIndex((1, 3))
julia> tri = grid.cells[cn]
Triangle((1, 2, 6))
julia> i, j, k = tri.nodes;
julia> ((i, j), (j, k), (k, i))[lcn]
(6, 1)
julia> grid.nodes[[6, 1]]
2-element Vector{Node{2, Float64}}:
Node{2, Float64}([-1.0, -0.5])
Node{2, Float64}([-1.0, -1.0])
shows that the associated triangle is (1, 2, 6)
and the edge (6, 1)
belongs to the left boundary. I'll try and write a simple example with this understanding.
Thanks for the input @DanielVandH ! Let me comment on this, trying to help you understand decisions made on our side.
Currently we essentially assume that we can query the following information from the grid efficiently:
grid.cells
grid.nodes
grid.cellsets
and so onFor subdomain information we use ordered grids to utilize caches as efficient as possible. If you use a Set of indices, then the iteration order does not follow the natural iteration order of your cells. In benchmarks we have shown that the iteration over the OrderedSet gives you quite a bit of speed up (4x and more, depending on how hard you try to optimize cache utilization).
For the indices itself we have the cells and nodes, which are just integers (as they index into grid.cells and grid.nodes directly) and we have furthermore indices for subsets of cells, e.g. facets. Since our grid does not materialize facets and ridges explicitly, we have for now settled on this trick to still be able to index into these entities.
If there are further questions do not hesitate to ask.
i, j, k = tri.nodes
((i, j), (j, k), (k, i))[lcn]
We have the function
Ferrite.facets(tri)[lcn]
which does this operation btw. As long as you only have linear triangles, this will return all nodes on the facet, as this corresponds to the vertices of the triangle. For 2nd order triangles, there is a difference between nodes (all points defining the triangle) and the vertices (only the corners of the triangles).
Thanks for your comments. Here are some functions that are useful for this. I might be using some internals of DelaunayTriangulation here, I'd have to double check. If I am, I'll make them public.
First, getting the cells and nodes is very easy. For get_cells
, I also a return a mapping that takes a triangle to its position in cells
, since my methods for getting the facetsets
would be pretty slow otherwise.
const DT = DelaunayTriangulation
function get_cells(tri::Triangulation)
cells = Vector{Triangle}(undef, num_solid_triangles(tri))
cell_map = Dict{NTuple{3,Int},Int}()
for (cell_number, triangle) in enumerate(each_solid_triangle(tri))
ijk = triangle_vertices(triangle)
uvw = DT.sort_triangle(ijk) # It's useful to know the order of the vertices
cells[cell_number] = Triangle(uvw)
cell_map[uvw] = cell_number
end
return cells, cell_map
end
function get_nodes(tri::Triangulation)
nodes = Vector{Node{2,Float64}}(undef, num_solid_vertices(tri))
for (node_number, vertex) in enumerate(each_point_index(tri)) # was originally each_solid_vertex, but that breaks if there are any points not in the triangulation.
p = get_point(tri, vertex)
x, y = getxy(p)
nodes[node_number] = Node((x, y))
end
return nodes
end
For getting the boundary information, there are actually several ways to do this inside DelaunayTriangulation.jl.
each_boundary_edge
. This is unordered.each_ghost_vertex
together with get_adjacent2vertex
. This is also unordered, but each of the facet sets can be built one at a time. get_boundary_nodes
. This requires a bit of care to get correct.I think the third is probably fastest to work with based on your comment @termi-official. Just to show what could be done, here are all three approaches. The third is the most involved but that's just how it is.
using OrderedCollections
function to_facetindex(cell_map, ijk, e)
uvw = DT.sort_triangle(ijk)
cell_number = cell_map[uvw]
u, v, w = triangle_vertices(uvw)
i, j = edge_vertices(e)
(i, j) == (u, v) && return FacetIndex(cell_number, 1)
(i, j) == (v, w) && return FacetIndex(cell_number, 2)
return FacetIndex(cell_number, 3) # assume (i, j) == (w, u)
end
function get_facetsets_1(tri::Triangulation, cell_map)
# Using each_boundary_edge (same as keys(get_boundary_edge_map(tri)))
facetsets = Dict{String,OrderedSet{FacetIndex}}()
ghost_map = Dict(i => string(i) for i in each_ghost_vertex(tri)) # trying to avoid allocating string(i) every single time
for e in DT.each_boundary_edge(tri)
u, v = edge_vertices(e)
w = get_adjacent(tri, e)
g = get_adjacent(tri, v, u) # ghost vertex
set = get!(OrderedSet{FacetIndex}, facetsets, ghost_map[g])
push!(set, to_facetindex(cell_map, (u, v, w), e))
end
return facetsets
end
function get_facetsets_2(tri::Triangulation, cell_map)
# Using each_ghost_vertex together with get_adjacent2vertex
facetsets = Dict{String,OrderedSet{FacetIndex}}()
for g in each_ghost_vertex(tri)
edges = get_adjacent2vertex(tri, g)
set = OrderedSet{FacetIndex}()
facetsets[string(g)] = set
for e′ in each_edge(edges)
e = DT.reverse_edge(e′) # e′ is the edge adjoining the ghost vertex
u, v = edge_vertices(e)
w = get_adjacent(tri, e)
push!(set, to_facetindex(cell_map, (u, v, w), e))
end
end
return facetsets
end
function get_facetsets_3(tri::Triangulation, cell_map)
# Iterating over the boundary nodes in order
facetsets = Dict{String,OrderedSet{FacetIndex}}()
g = 0
nc = DT.num_curves(tri)
for curve_index in 1:nc
if nc == 1
bn = get_boundary_nodes(tri)
else
bn = get_boundary_nodes(tri, curve_index)
end
ns = DT.num_sections(bn)
for segment_index in 1:ns
g -= 1
set = OrderedSet{FacetIndex}()
facetsets[string(g)] = set
if nc == ns == 1
bnn = bn
else
bnn = get_boundary_nodes(bn, segment_index)
end
ne = num_boundary_edges(bnn)
for i in 1:ne
u = get_boundary_nodes(bnn, i)
v = get_boundary_nodes(bnn, i + 1)
w = get_adjacent(tri, u, v)
push!(set, to_facetindex(cell_map, (u, v, w), (u, v)))
end
end
end
return facetsets
end
A function to then generate a grid would be
function Ferrite.Grid(tri::Triangulation)
cells, cell_map = get_cells(tri)
nodes = get_nodes(tri)
facetsets = get_facetsets_3(tri, cell_map)
return Ferrite.Grid(cells, nodes; facetsets)
end
As an example, here I create a Grid
for an annulus.
function Annulus(R₁, R₂) # centred at the origin
@assert R₁ < R₂
inner = CircularArc((R₁, 0.0), (R₁, 0.0), (0.0, 0.0), positive = false)
outer = CircularArc((R₂, 0.0), (R₂, 0.0), (0.0, 0.0))
return [[[outer]], [[inner]]]
end
R₁, R₂ = 1.0, 2.0
tri = triangulate(NTuple{2, Float64}[], boundary_nodes = Annulus(R₁, R₂))
refine!(tri; max_area = 1e-3π * (R₂^2 - R₁^2))
julia> grid = Grid(tri)
Grid{2, Triangle, Float64} with 1614 Triangle cells and 890 nodes
julia> grid.facetsets
Dict{String, OrderedSet{FacetIndex}} with 2 entries:
"-2" => OrderedSet{FacetIndex}([FacetIndex((1473, 3)), FacetIndex((129, 2)), FacetIndex((257, 3)), FacetIndex((196, 2)), FacetIndex((1534, 3)), FacetIndex((1295, 3)), FacetIndex((671, 2)), FacetInde…
"-1" => OrderedSet{FacetIndex}([FacetIndex((686, 3)), FacetIndex((1469, 3)), FacetIndex((105, 2)), FacetIndex((331, 3)), FacetIndex((673, 2)), FacetIndex((252, 3)), FacetIndex((1463, 2)), FacetIndex…
What do you think?
I handle the boundary conditions a bit differently inside FiniteVolumeMethod.jl (https://github.com/SciML/FiniteVolumeMethod.jl/blob/9dc2612219b95eb49b4a81d8cace441f10b6ded4/src/conditions.jl#L452C1-L490C1), which also allows for internal constraints (useful for e.g. this example https://sciml.github.io/FiniteVolumeMethod.jl/stable/wyos/poissons_equation/#Using-the-Provided-Template). I don't know if you allow for internal constraints in this package-I imagine you do-, but it wouldn't be too hard for me to include that. Probably out of scope for this issue but just as an FYI I guess.
Might be something wrong in my code above, not really sure. I tried using it for https://ferrite-fem.github.io/Ferrite.jl/dev/tutorials/heat_equation/, using triangles instead of quadrilaterals but I get some singular matrix:
using Ferrite, SparseArrays, DelaunayTriangulation
tri = triangulate_rectangle(-1, 1, -1, 1, 20, 20, single_boundary=true)
grid = Grid(tri)
ip = Lagrange{RefTriangle,1}()
qr = QuadratureRule{RefTriangle}(2)
cellvalues = CellValues(qr, ip)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
K = allocate_matrix(dh)
ch = ConstraintHandler(dh)
∂Ω = getfacetset(grid, "-1")
dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0.0)
add!(ch, dbc)
close!(ch)
function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues)
n_basefuncs = getnbasefunctions(cellvalues)
# Reset to 0
fill!(Ke, 0)
fill!(fe, 0)
# Loop over quadrature points
for q_point in 1:getnquadpoints(cellvalues)
# Get the quadrature weight
dΩ = getdetJdV(cellvalues, q_point)
# Loop over test shape functions
for i in 1:n_basefuncs
δu = shape_value(cellvalues, q_point, i)
∇δu = shape_gradient(cellvalues, q_point, i)
# Add contribution to fe
fe[i] += δu * dΩ
# Loop over trial shape functions
for j in 1:n_basefuncs
∇u = shape_gradient(cellvalues, q_point, j)
# Add contribution to Ke
Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
end
end
end
return Ke, fe
end
function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler)
# Allocate the element stiffness matrix and element force vector
n_basefuncs = getnbasefunctions(cellvalues)
Ke = zeros(n_basefuncs, n_basefuncs)
fe = zeros(n_basefuncs)
# Allocate global force vector f
f = zeros(ndofs(dh))
# Create an assembler
assembler = start_assemble(K, f)
# Loop over all cels
for cell in CellIterator(dh)
# Reinitialize cellvalues for this cell
reinit!(cellvalues, cell)
# Compute element contribution
assemble_element!(Ke, fe, cellvalues)
# Assemble Ke and fe into K and f
assemble!(assembler, celldofs(cell), Ke, fe)
end
return K, f
end
K, f = assemble_global(cellvalues, K, dh)
ERROR: ArgumentError: det(J) is not positive: det(J) = -1.35180055401662
Stacktrace:
[1] throw_detJ_not_pos(detJ::Float64)
@ Ferrite C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\common_values.jl:5
[2] _update_detJdV!
@ C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\CellValues.jl:96 [inlined]
[3] reinit!(cv::CellValues{…}, cell::Nothing, x::Vector{…})
@ Ferrite C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\CellValues.jl:119
[4] reinit!
@ C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\CellValues.jl:102 [inlined]
[5] reinit!
@ C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\iterators.jl:92 [inlined]
[6] assemble_global(cellvalues::CellValues{…}, K::SparseMatrixCSC{…}, dh::DofHandler{…})
@ Main .\Untitled-1:186
[7] top-level scope
@ Untitled-1:195
Some type information was truncated. Use `show(err)` to see complete types.
Will try look some more into it, but just posting for now in case there's anything obvious I'm missing.
That is probably because Ferrite requires the node order to be counter clockwise. (Just skimming through, I guess the sorting of the node ids (DT.sort_triangl
) could be the problem?)
sort_triangle
preserves the order of the triangle (counter-clockwise), it just rotates it to put the minimum vertex first. I think the orientation is consistent throughout but I'll check.
for T in grid.cells
ijk = T.nodes
DT.is_positively_oriented(DT.triangle_orientation(tri, ijk)) || error("Triangle is negatively oriented")
end
shows it's fine.
Nevermind. It's because I fixed an issue in my post above but somehow didn't do the same in my actual script 🤦 .
Here's a complete example. It's https://ferrite-fem.github.io/Ferrite.jl/dev/tutorials/heat_equation/ except on an annulus. The boundary conditions are zero on the inner boundary, and $x^2y^2$ on the outer boundary. I use tricontourf
to plot the solution.
using Ferrite, SparseArrays, CairoMakie, DelaunayTriangulation, OrderedCollections
const DT = DelaunayTriangulation
function get_cells(tri::Triangulation)
cells = Vector{Triangle}(undef, num_solid_triangles(tri))
cell_map = Dict{NTuple{3,Int},Int}()
for (cell_number, triangle) in enumerate(each_solid_triangle(tri))
ijk = triangle_vertices(triangle)
uvw = DT.sort_triangle(ijk) # It's useful to know the order of the vertices
cells[cell_number] = Triangle(uvw)
cell_map[uvw] = cell_number
end
return cells, cell_map
end
function get_nodes(tri::Triangulation)
nodes = Vector{Node{2,Float64}}(undef, DT.num_points(tri))
for vertex in DT.each_point_index(tri)
p = get_point(tri, vertex)
x, y = getxy(p)
nodes[vertex] = Node((x, y))
end
return nodes
end
function to_facetindex(cell_map, ijk, e)
uvw = DT.sort_triangle(ijk)
cell_number = cell_map[uvw]
u, v, w = triangle_vertices(uvw)
i, j = edge_vertices(e)
(i, j) == (u, v) && return FacetIndex(cell_number, 1)
(i, j) == (v, w) && return FacetIndex(cell_number, 2)
return FacetIndex(cell_number, 3) # assume (i, j) == (w, u)
end
function get_facetsets(tri::Triangulation, cell_map)
# Iterating over the boundary nodes in order
facetsets = Dict{String,OrderedSet{FacetIndex}}()
g = 0
nc = DT.num_curves(tri)
for curve_index in 1:nc
if nc == 1
bn = get_boundary_nodes(tri)
else
bn = get_boundary_nodes(tri, curve_index)
end
ns = DT.num_sections(bn)
for segment_index in 1:ns
g -= 1
set = OrderedSet{FacetIndex}()
facetsets[string(g)] = set
if nc == ns == 1
bnn = bn
else
bnn = get_boundary_nodes(bn, segment_index)
end
ne = num_boundary_edges(bnn)
for i in 1:ne
u = get_boundary_nodes(bnn, i)
v = get_boundary_nodes(bnn, i + 1)
w = get_adjacent(tri, u, v)
push!(set, to_facetindex(cell_map, (u, v, w), (u, v)))
end
end
end
return facetsets
end
function Ferrite.Grid(tri::Triangulation)
cells, cell_map = get_cells(tri)
nodes = get_nodes(tri)
facetsets = get_facetsets(tri, cell_map)
return Ferrite.Grid(cells, nodes; facetsets)
end
function Annulus(R₁, R₂) # centred at the origin
@assert R₁ < R₂
inner = CircularArc((R₁, 0.0), (R₁, 0.0), (0.0, 0.0), positive=false)
outer = CircularArc((R₂, 0.0), (R₂, 0.0), (0.0, 0.0))
return [[[outer]], [[inner]]]
end
R₁, R₂ = 1.0, 2.0
tri = triangulate(NTuple{2,Float64}[], boundary_nodes=Annulus(R₁, R₂))
refine!(tri; max_area=1e-3π * (R₂^2 - R₁^2))
grid = Grid(tri)
ip = Lagrange{RefTriangle,1}()
qr = QuadratureRule{RefTriangle}(2)
cellvalues = CellValues(qr, ip)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
K = allocate_matrix(dh)
ch = ConstraintHandler(dh)
dbc = Dirichlet(:u, getfacetset(grid, "-1"), (x, t) -> x[1]^2 * x[2]^2)
add!(ch, dbc)
nbc = Dirichlet(:u, getfacetset(grid, "-2"), (x, t) -> 0.0)
add!(ch, nbc)
close!(ch)
function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues)
n_basefuncs = getnbasefunctions(cellvalues)
# Reset to 0
fill!(Ke, 0)
fill!(fe, 0)
# Loop over quadrature points
for q_point in 1:getnquadpoints(cellvalues)
# Get the quadrature weight
dΩ = getdetJdV(cellvalues, q_point)
# Loop over test shape functions
for i in 1:n_basefuncs
δu = shape_value(cellvalues, q_point, i)
∇δu = shape_gradient(cellvalues, q_point, i)
# Add contribution to fe
fe[i] += δu * dΩ
# Loop over trial shape functions
for j in 1:n_basefuncs
∇u = shape_gradient(cellvalues, q_point, j)
# Add contribution to Ke
Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
end
end
end
return Ke, fe
end
function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler)
# Allocate the element stiffness matrix and element force vector
n_basefuncs = getnbasefunctions(cellvalues)
Ke = zeros(n_basefuncs, n_basefuncs)
fe = zeros(n_basefuncs)
# Allocate global force vector f
f = zeros(ndofs(dh))
# Create an assembler
assembler = start_assemble(K, f)
# Loop over all cels
for cell in CellIterator(dh)
# Reinitialize cellvalues for this cell
reinit!(cellvalues, cell)
# Compute element contribution
assemble_element!(Ke, fe, cellvalues)
# Assemble Ke and fe into K and f
assemble!(assembler, celldofs(cell), Ke, fe)
end
return K, f
end
K, f = assemble_global(cellvalues, K, dh)
apply!(K, f, ch)
u = K \ f
points = [p.x for p in getnodes(grid)]
projector = L2Projector(ip, grid)
ph = PointEvalHandler(grid, points)
sol = evaluate_at_points(ph, projector, u)
tricontourf(tri, sol)
Let me know if there is anything else you think could be useful to demonstrate.
That looks very nice! If you want to, it would be great to open a PR adding this as a how-to IMO. Then it is also easy to comment / test it!
Regarding the postprocessing, I'd suggest
- points = [p.x for p in getnodes(grid)]
- projector = L2Projector(ip, grid)
- ph = PointEvalHandler(grid, points)
- sol = evaluate_at_points(ph, projector, u)
+ sol = evaluate_at_grid_nodes(dh, u, :u)
(The fact that evaluate_at_points(ph, projector, u)
works is a bit unfortunate in this case actually, since u
then should be the result of a projection and be numbered according to the projector's DofHandler, but I don't know how we could detect this. In this case, it will give the correct results since the dof handlers are identical.)
Thanks for the tip with evaluate_at_grid_nodes
! I figured there was an easier way but that was the best I could find in the docs.
I'll write up the example properly with Literate and make a PR soon
It would be nice to have an example where https://github.com/JuliaGeometry/DelaunayTriangulation.jl is used for the grid generation. Perhaps we can spice up one of the existing examples that uses a hypercube grid to something more fun.
The author of DelaunayTriangulation.jl (@DanielVandH) have offered to help out with any questions that may arise.