Refinement of getcoordinates and addition of getvolume

[zkk960317]

There is a lack of functions to get the volume of in Ferrite. Here I have given a solution. It needs to be checked by the developer and considered for merging. The following is the functions defined outside of module Ferrite. If put these to grid.jl, the prefix Ferrite. should be removed.

function Ferrite.getcoordinates(grid::Ferrite.AbstractGrid, nodes::NTuple{N,Int64}) where {N}
    dim = Ferrite.getdim(grid)
    T = Ferrite.get_coordinate_eltype(grid)
    n = lastindex(nodes)
    x = Vector{Vec{dim,T}}(undef, n)
    @inbounds for i in 1:n
        x[i] = Ferrite.getcoordinates(getnodes(grid, nodes[i]))
    end
    return x
end

## The  following old function in Ferrite 0.3.14 is to get the coordinates of all nodes of a cell.
## function getcoordinates(grid::AbstractGrid, cell::CellIndex) = (grid, grid.cells[cell.idx[1]].nodes)
## get the coordinates of all vertices of cell
Ferrite.getcoordinates(grid::Ferrite.AbstractGrid, cell::CellIndex) = Ferrite.getcoordinates(grid, Ferrite.vertices(grid.cells[cell.idx[1]]))
Ferrite.getcoordinates(grid::Ferrite.AbstractGrid, face::FaceIndex) = getcoordinates(grid, Ferrite.faces(grid.cells[face[1]])[face[2]])
Ferrite.getcoordinates(grid::Ferrite.AbstractGrid, edge::EdgeIndex) = getcoordinates(grid, Ferrite.edges(grid.cells[edge[1]])[edge[2]])
Ferrite.getcoordinates(grid::Ferrite.AbstractGrid, vertex::VertexIndex) = getcoordinates(grid, (Ferrite.vertices(grid.cells[vertex[1]])[vertex[2]],))


# The `getvolume` function of a 3d / 2d / 1d cell is defined to get its volume / area / length, respectively. 
# Specially, the volume of 3d planar cell Quadrilateral3D is zero, FaceIndex(x,1) should be used as the second parameter to get the area of Quadrilateral3D cell x. 
# Similarly, the volume of Line3D and Line2D is zero, the correct parameter is EdgeIndex(x,1) and FaceIndex(x,1).
getvolume(grid::Ferrite.AbstractGrid{3}, cell::CellIndex) = _getvolume(Ferrite.getcoordinates(grid, cell))
getvolume(grid::Ferrite.AbstractGrid{3}, face::FaceIndex) = _getarea(Ferrite.getcoordinates(grid, face))
getvolume(grid::Ferrite.AbstractGrid{3}, edge::EdgeIndex) = _getlength(Ferrite.getcoordinates(grid, edge))
getvolume(grid::Ferrite.AbstractGrid{2}, cell::CellIndex) = _getarea(Ferrite.getcoordinates(grid, cell))
getvolume(grid::Ferrite.AbstractGrid{2}, face::FaceIndex) = _getlength(Ferrite.getcoordinates(grid, face))
getvolume(grid::Ferrite.AbstractGrid{1}, cell::CellIndex) = _getlength(Ferrite.getcoordinates(grid, cell))
#alias
getvolume(grid::Ferrite.AbstractGrid{3}, cell::Int) = getvolume(grid::Ferrite.AbstractGrid{3}, CellIndex(cell))

_getlength(coords::Vector{Vec{dim,T}}) where {dim,T} = norm(coords[1] - coords[2])

# Assuming that the faces of any 3D cell are planar.
# 2-order cell are calculated as 1-order cell (such as QuadraticTetrahedron -> Tetrahedron), ignoring the nodes in the edges.
# Assuming that the coordinate points enclose a convex area in counterclockwise order.
function _getarea(coords::Vector{Vec{dim,T}}) where {dim,T}
    n = lastindex(coords)
    n < 3 && return 0.0
    vec_coords = [coords[i] - coords[1] for i in 2:n]
    area = 0.0
    for i in 1:n-2
        area += norm(cross(vec_coords[i], vec_coords[i+1])) / 2
    end
    return area
end

function _getvolume(coords::Vector{Vec{dim,T}}) where {dim,T}
    @inline volu_tet(coords) = dot(cross(coords[2] - coords[1], coords[3] - coords[1]), coords[4] - coords[1]) / 6
    n = lastindex(coords)
    n < 4 && return 0.0
    volu = 0.0
    if n == 4 # Tetrahedron
        volu = volu_tet(coords)
    elseif n == 8 # Hexahedron, divided into 5 Tetrahedron.
        volu = volu_tet(coords[[2, 1, 6, 3]]) + volu_tet(coords[[7, 3, 6, 8]]) +
               volu_tet(coords[[4, 1, 3, 8]]) + volu_tet(coords[[5, 1, 8, 6]]) +
               volu_tet(coords[[6, 1, 8, 3]])
    else
        error("The volume of $(dim)D cell with $n vertices is not supported")
    end
    return volu
end

The following is the test result:

# 2d test
nodes_2d = [Node((0.0, 0.0)), Node((1.0, 0.0)), Node((1.0, 1.0)), Node((0.0, 1.0))]
cells_2d = [Line2D((1, 2)), Triangle((1, 2, 3)), Quadrilateral((1, 2, 3, 4))]
grid = Grid(cells_2d, nodes_2d)
#Line2D
getvolume(grid, CellIndex(1)) # 0.0
getvolume(grid, FaceIndex(1, 1)) # 1.0
getvolume(grid, EdgeIndex(1, 1)) # ERROR!
#Triangle
getvolume(grid, CellIndex(2)) # 0.5
getvolume(grid, FaceIndex(2, 1)) # 1.0
getvolume(grid, EdgeIndex(2, 1)) # ERROR! 
#Quadrilateral
getvolume(grid, CellIndex(3)) # 1.0
getvolume(grid, FaceIndex(3, 1)) # 1.0
getvolume(grid, EdgeIndex(3, 1)) # ERROR! 

# 3d test
nodes_3d = [Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)), Node((1.0, 1.0, 0.0)), Node((0.0, 1.0, 0.0)),
    Node((0.0, 0.0, 1.0)), Node((1.0, 0.0, 1.0)), Node((1.0, 1.0, 1.0)), Node((0.0, 1.0, 1.0))]
cells_3d = [Line3D((1, 2)), Quadrilateral3D((1, 2, 3, 4)), Tetrahedron((1, 2, 3, 5)), Hexahedron((1, 2, 3, 4, 5, 6, 7, 8))]
grid = Grid(cells_3d, nodes_3d)
#Line3D
getvolume(grid, CellIndex(1)) # 0.0
getvolume(grid, FaceIndex(1, 1)) # ERROR!
getvolume(grid, EdgeIndex(1, 1)) # 1.0
#Quadrilateral3D
getvolume(grid, CellIndex(2)) # 0.0
getvolume(grid, FaceIndex(2, 1)) # 1.0
getvolume(grid, EdgeIndex(2, 1)) # 1.0
#Tetrahedron
getvolume(grid, CellIndex(3)) # 0.16666
getvolume(grid, FaceIndex(3, 3)) # 0.7071
getvolume(grid, EdgeIndex(3, 1)) # 1.0
#Hexahedron
getvolume(grid, CellIndex(4)) # 1.0
getvolume(grid, FaceIndex(4, 1)) # 1.0
getvolume(grid, EdgeIndex(4, 1)) # 1.0

[termi-official]

What exactly is the question or issue?

[zkk960317]

@termi-official The issue is lack of functions to get volume of a cell in Ferrite, which is important. Here I give a solution. Need developers to check and consider merging.

[KnutAM]

Hi and thanks for your suggestions.
Note that you can get the volume of a cell by using the CellValues, e.g.

using Ferrite
# General setup
grid = generate_grid(Quadrilateral, (2,3))
qr = QuadratureRule{RefQuadrilateral}(1)
ip = geometric_interpolation(getcelltype(grid))
cv = CellValues(qr, ip, ip)
# Get volume of cell nr 2
cellnr = 2
reinit!(cv, getcoordinates(grid, cellnr))
volume = sum(getdetJdV(cv, i) for i in 1:getnquadpoints(cv))
# Or in this case, we only have one quadrature points, so the following works as well
volume2 = getdetJdV(cv, 1)

The possibility of getting coordinates for specific entities, such as faces, is a good idea. But in case of higher order cells (with nonlinear geometric interpolations), the above won't work. See some related discussion here #702.