Document grid.jl

[mforets]

In the process of getting to know the library, i'd like to contribute with the docs for Grid.

Here is a list (freel free to edit):

• [x] Describe grid format. Node, Cell and Grid.
• [x] Grid utility functions.
• [ ] Grid transformations.
• [ ] Types of cells.
• [x] How to add sets.
• [ ] Describe mixed grids.

[KristofferC]

Describing the way a global face is defined ((global_cell_index, local_face_index)) would also be good. The local_face_index is defined by e.g.

https://github.com/KristofferC/JuAFEM.jl/blob/8224282ab4d67cb523ef342e4a6ceb1716764ada/src/interpolations.jl#L154

To give a description of it, say that we have a "surface" (line) on a 2D mesh, represented by the global nodes 6, 10. And let's say that we have an element with number 3 and with the global vs local nodes:

15----10     4-----3
|     |      |     |
| (3) |      |     |
|     |      |     |
3-----6      1-----2

global        local

Since the global nodes (6, 10) are at the side with the local nodes (2, 3) and the local surface (2,3) has index 2 in the function

https://github.com/KristofferC/JuAFEM.jl/blob/8224282ab4d67cb523ef342e4a6ceb1716764ada/src/interpolations.jl#L154

the surface with global nodes 6 and 10 would be stored as (3, 2) (global cell 3 with local face 2).

[KristofferC]

I updated the comment above with a bit more details.

[KristofferC]

Here is a an algorithm for computing the JuAFEM datastructure of the faces when you have a 2D mesh and a lower dimensional "face mesh". I believe it should be efficient enough.

The example I am using it on is:

using JuAFEM

# elements from example 2d mesh above
elements = [
    (1, 2, 5),
    (2, 3, 6),
    (1, 5, 4),
    (2, 6, 5),
    (4, 5, 8),
    (5, 6, 9),
    (4, 8, 7),
    (5, 9, 8)
]

# edges from example face mesh above
edges = [
    (3, 6),
    (6, 9)
]

# The "interpolation" we are using for the mesh, here 2D linear triangles
interpolation = Lagrange{2, RefTetrahedron, 1}()

# algorithm to compute the faces
function compute_faceset(elements, edges, ip::Interpolation{dim}) where {dim}
    local_faces = JuAFEM.faces(ip)
    nodes_per_face = length(local_faces[1])
    d = Dict{NTuple{nodes_per_face, Int}, Tuple{Int, Int}}()
    for (e, element) in enumerate(elements) # e is global element number
        for (f, face) in enumerate(local_faces) # f is local face number
            # store the global nodes for the particular element, local face combination
            d[ntuple(i-> element[face[i]], nodes_per_face)] = (e, f)
        end
    end

    faces = Set{Tuple{Int, Int}}()
    for edge in edges
        # lookup the element, local face combination for this edge
        push!(faces, d[edge])
    end

    return faces
end

Running it:

julia> compute_cell_local_edge(elements, edges, interpolation)
Set{Tuple{Int64,Int64}} with 2 elements:
  (6, 2)
  (2, 2)

And indeed, the two faces are for the local face number 2 on cells 2 and 6.

Another example:

julia> edges_2 = [
           (9, 8),
           (8, 7),
           (7, 4),
           (4, 1),
       ];

julia> compute_cell_local_edge(elements, edges_2, interpolation)
Set{Tuple{Int64,Int64}} with 4 elements:
  (8, 2)
  (7, 2)
  (3, 3)
  (7, 3)

[jorgepz]

Perfect! I've extended the example to test it in a 3D mesh. And it works very well! I considered a cube with 8 nodes and 6 tetrahedrons within that cube

The code of the example is

using JuAFEM

# elements from example 2d mesh above
elements2DMesh = [
    (1, 2, 5),
    (2, 3, 6),
    (1, 5, 4),
    (2, 6, 5),
    (4, 5, 8),
    (5, 6, 9),
    (4, 8, 7),
    (5, 9, 8)
]

# edges from example face mesh above
edges2DMesh = [
    (3, 6),
    (6, 9)
]

# The "interpolation" we are using for the mesh, here 2D linear triangles
interpolation2DMesh = Lagrange{2, RefTetrahedron, 1}()

# elements from example 3d mesh above
elements3DMesh = [
  (1, 4, 2, 6), 
  (6, 2, 3, 4),
  (4, 3, 6, 7),
  (4, 1, 5, 6),
  (4, 6, 5, 8),
  (4, 7, 6, 8)
]

# edges from example face mesh above
edges3DMeshA = [
    (7, 6, 8),
    (6, 5, 8)
]

edges3DMeshB = [
    (6, 7, 8),
    (6, 5, 8)
]

# The "interpolation" we are using for the mesh, here 3D linear tetrahedrons
interpolation3DMesh = Lagrange{3, RefTetrahedron, 1}()

# algorithm to compute the faces
function compute_faceset(elements, edges, ip::Interpolation{dim}) where {dim}
    local_faces = JuAFEM.faces(ip)
    nodes_per_face = length(local_faces[1])
    d = Dict{NTuple{nodes_per_face, Int}, Tuple{Int, Int}}()
    for (e, element) in enumerate(elements) # e is global element number
        for (f, face) in enumerate(local_faces) # f is local face number
            # store the global nodes for the particular element, local face combination
            d[ntuple(i-> element[face[i]], nodes_per_face)] = (e, f)
        end
    end

    faces = Set{Tuple{Int, Int}}()
    for edge in edges
        # lookup the element, local face combination for this edge
        push!(faces, d[edge])
    end

    return faces
end

there are two edges. In the case where the orientation of the surface triangles is the same as the nodes in the tetrahedron it works fine

julia> compute_faceset(elements3DMesh, edges3DMeshA, interpolation3DMesh)

Set(Tuple{Int64,Int64}[(6, 3), (5, 3)])

if not, it does not work since it does not find an exact match for the triangle.

julia> compute_faceset(elements3DMesh, edges3DMeshB, interpolation3DMesh)
ERROR: KeyError: key (6, 7, 8) not found
Stacktrace:
 [1] getindex(::Dict{Tuple{Int64,Int64,Int64},Tuple{Int64,Int64}}, ::Tuple{Int64,Int64,Int64}) at ./dict.jl:478
 [2] compute_faceset(::Array{NTuple{4,Int64},1}, ::Array{Tuple{Int64,Int64,Int64},1}, ::Lagrange{3,RefTetrahedron,1}) at /home/jor/testFaces.jl:70
 [3] top-level scope at none:0

The orientation could be of interest if for instance pressure loads are given in a local coordinate system with the orientation given in the mesh. It might be of interest to store the orientation as well.