BiVeCo + VTK reader (#93)

Project.toml

@@ -15,6 +15,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
+ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"

docs/Manifest.toml

@@ -478,8 +478,8 @@ uuid = "29a986be-02c6-4525-aec4-84b980013641"
 version = "2.0.4"
 
 [[deps.Ferrite]]
-deps = ["EnumX", "LinearAlgebra", "NearestNeighbors", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
-git-tree-sha1 = "251e85da52ac1ae1944265af9c42629df72c61e0"
+deps = ["EnumX", "ForwardDiff", "LinearAlgebra", "NearestNeighbors", "OrderedCollections", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
+git-tree-sha1 = "e608f28615d999ae7be3172e5af186236030a0d0"
 repo-rev = "master"
 repo-url = "https://github.com/Ferrite-FEM/Ferrite.jl.git"
 uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
@@ -1366,6 +1366,12 @@ git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111"
 uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
 version = "1.5.3"
 
+[[deps.ReadVTK]]
+deps = ["Base64", "CodecZlib", "Downloads", "LightXML", "Reexport", "VTKBase"]
+git-tree-sha1 = "f8a48e99ca616b46ad62356257dd21b9bb522024"
+uuid = "dc215faf-f008-4882-a9f7-a79a826fadc3"
+version = "0.2.0"
+
 [[deps.RecipesBase]]
 deps = ["PrecompileTools"]
 git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
@@ -1772,7 +1778,7 @@ uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
 version = "0.5.2"
 
 [[deps.Thunderbolt]]
-deps = ["BlockArrays", "CommonSolve", "Ferrite", "FerriteGmsh", "JLD2", "Krylov", "LinearAlgebra", "ModelingToolkit", "OrderedCollections", "Polyester", "Reexport", "SparseArrays", "SparseMatricesCSR", "StaticArrays", "SymbolicIndexingInterface", "Tensors", "TimerOutputs", "UnPack"]
+deps = ["BlockArrays", "CommonSolve", "DataStructures", "DiffEqBase", "Ferrite", "FerriteGmsh", "JLD2", "Krylov", "LinearAlgebra", "ModelingToolkit", "OrderedCollections", "Polyester", "ReadVTK", "Reexport", "SparseArrays", "SparseMatricesCSR", "StaticArrays", "SymbolicIndexingInterface", "Tensors", "TimerOutputs", "UnPack", "Unrolled", "WriteVTK"]
 path = ".."
 uuid = "909927c2-98d5-4a67-bba9-79f03a9ad49b"
 version = "0.0.1-DEV"
@@ -1848,6 +1854,12 @@ git-tree-sha1 = "25008b734a03736c41e2a7dc314ecb95bd6bbdb0"
 uuid = "a7c27f48-0311-42f6-a7f8-2c11e75eb415"
 version = "0.1.6"
 
+[[deps.Unrolled]]
+deps = ["MacroTools"]
+git-tree-sha1 = "6cc9d682755680e0f0be87c56392b7651efc2c7b"
+uuid = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
+version = "0.1.5"
+
 [[deps.VTKBase]]
 git-tree-sha1 = "c2d0db3ef09f1942d08ea455a9e252594be5f3b6"
 uuid = "4004b06d-e244-455f-a6ce-a5f9919cc534"

src/Thunderbolt.jl

@@ -11,7 +11,8 @@ using OrderedCollections
 using BlockArrays, SparseArrays, StaticArrays
 
 using JLD2
-using WriteVTK
+import WriteVTK
+import ReadVTK
 
 # This is a standalone module which will be a custom package in the future
 include("solver/operator_splitting.jl")

src/io.jl

@@ -208,3 +208,82 @@ end
 
 finalize_timestep!(io::JLD2Writer, t) = nothing
 finalize!(io::JLD2Writer) = nothing
+
+function reorder_nodal!(dh::DofHandler)
+    @assert length(dh.field_names) == 1 "Just single field possible."
+    grid = Ferrite.get_grid(dh)
+    for sdh in dh.subdofhandlers
+        firstcellidx = first(sdh.cellset)
+        celltype = typeof(getcells(grid,firstcellidx))
+        @assert sdh.field_interpolations[1] == Ferrite.geometric_interpolation(celltype)
+        for i ∈ 1:getncells(grid)
+            dh.cell_dofs[dh.cell_dofs_offset[i]:(dh.cell_dofs_offset[i]+Ferrite.ndofs_per_cell(dh, i)-1)] .= getcells(grid,i).nodes
+        end
+    end
+end
+
+function to_ferrite_elements(cells_vtk::Vector{WriteVTK.MeshCell{WriteVTK.VTKCellType, Vector{Int64}}})
+    celltype = if cells_vtk[1].ctype == WriteVTK.VTKCellTypes.VTK_TETRA
+        Tetrahedron
+    elseif cells_vtk[1].ctype == WriteVTK.VTKCellTypes.VTK_HEXAHEDRON
+        Hexahedron
+    else 
+        @error "Unknown cell type" cells_vtk[1].ctype
+    end
+    cells_ferrite = Vector{celltype}(undef, length(cells_vtk))
+    for (i,cell_vtk) in enumerate(cells_vtk)
+        cells_ferrite[i] = if cells_vtk[1].ctype == WriteVTK.VTKCellTypes.VTK_TETRA
+            Tetrahedron(ntuple(i->cell_vtk.connectivity[i],4))
+        elseif cells_vtk[1].ctype == WriteVTK.VTKCellTypes.VTK_HEXAHEDRON
+            Hexahedron(ntuple(i->cell_vtk.connectivity[i],8))
+        end
+    end
+    return cells_ferrite
+end
+
+function read_vtk_cobivec(filename::String, transmural_id::String, apicobasal_id::String, radial_id::String, transventricular_id::String)
+    vtk = ReadVTK.VTKFile(filename)
+    points_vtk = get_points(vtk)
+    points_ferrite = [Vec{3}(point) for point in eachcol(points_vtk)]
+    cells_vtk = to_meshcells(get_cells(vtk))
+    cells_ferrite = to_ferrite_elements(cells_vtk)
+
+    grid = Grid(cells_ferrite, Node.(points_ferrite))
+
+    gip =  Ferrite.geometric_interpolation(typeof(grid.cells[1]))
+
+    dh = DofHandler(grid)
+    add!(dh, :u, gip)
+    close!(dh)
+    reorder_nodal!(dh)
+
+    all_data = get_point_data(vtk)
+    u_transmural = get_data(all_data[transmural_id])
+    u_apicobasal = get_data(all_data[apicobasal_id])
+    u_radial = get_data(all_data[radial_id])
+    u_transventricular = get_data(all_data[transventricular_id])
+
+    epicardium = OrderedSet{FacetIndex}()
+    endocardium = OrderedSet{FacetIndex}()
+    for (cellidx,cell) ∈ enumerate(grid.cells)
+        for (faceidx,facenodes) in enumerate(Ferrite.faces(cell))
+            # facenodes = cell.nodes[collect(face)]
+            if all(u_transmural[collect(facenodes)] .> 1.0-1e-6)
+                push!(endocardium, FacetIndex(cellidx, faceidx))
+            end
+            if all(u_transmural[collect(facenodes)] .< 1e-6)
+                push!(epicardium, FacetIndex(cellidx, faceidx))
+            end
+        end
+    end
+    Ferrite.addfacetset!(grid, "Epicardium", epicardium)
+    Ferrite.addfacetset!(grid, "Endocardium", endocardium)
+
+    return BiVCoordinateSystem(
+        dh,
+        collect(u_transmural),
+        collect(u_apicobasal),
+        collect(u_radial),
+        collect(u_transventricular)
+    )
+end

src/mesh/coordinate_systems.jl

@@ -253,3 +253,43 @@ function vtk_coordinate_system(vtk, cs::LVCoordinateSystem)
     vtk_point_data(vtk, cs.dh, cs.u_apicobasal, "apicobasal_")
     vtk_point_data(vtk, cs.dh, cs.u_transmural, "transmural_")
 end
+
+"""
+    BiVCoordinateSystem(dh, u_transmural, u_apicobasal, u_rotational, u_transventricular)
+
+Universal ventricular coordinate, containing the transmural, apicobasal, circumferential 
+and transventricular coordinates.
+"""
+struct BiVCoordinateSystem{DH <: Ferrite.AbstractDofHandler}
+    dh::DH
+    u_transmural::Vector{Float32}
+    u_apicobasal::Vector{Float32}
+    u_rotational::Vector{Float32}
+    u_transventricular::Vector{Float32}
+end
+
+
+"""
+BiVCoordinate{T}
+
+Biventricular universal coordinate, containing
+    * transmural
+    * apicobasal
+    * rotational
+    * transventricular
+"""
+struct BiVCoordinate{T}
+    transmural::T
+    apicaobasal::T
+    rotational::T
+    transventricular::T
+end
+
+getcoordinateinterpolation(cs::BiVCoordinateSystem, cell::Ferrite.AbstractCell) = Ferrite.getfieldinterpolation(cs.dh, (1,1))
+
+function vtk_coordinate_system(vtk, cs::BiVCoordinateSystem)
+    vtk_point_data(vtk, bivcs.dh, bivcs.u_transmural, "_transmural")
+    vtk_point_data(vtk, bivcs.dh, bivcs.u_apicobasal, "_apicobasal")
+    vtk_point_data(vtk, bivcs.dh, bivcs.u_rotational, "_rotational")
+    vtk_point_data(vtk, bivcs.dh, bivcs.u_transventricular, "_transventricular")
+end

src/modeling/core/coefficients.jl

@@ -111,6 +111,24 @@ function _evaluate_coefficient(coeff::CoordinateSystemCoefficient{<:LVCoordinate
     return LVCoordinate(x[1], x[2], x[3])
 end
 
+function evaluate_coefficient(coeff::CoordinateSystemCoefficient{<:BiVCoordinateSystem}, cell_cache, qp::QuadraturePoint{ref_shape,T}, t) where {ref_shape,T}
+    ip = getcoordinateinterpolation(coeff.cs, getcells(cell_cache.grid, cellid(cell_cache)))
+    dofs = celldofsview(coeff.cs.dh, cellid(cell_cache))
+    return _evaluate_coefficient(coeff, ip, dofs, qp, t)
+end
+
+function _evaluate_coefficient(coeff::CoordinateSystemCoefficient{<:BiVCoordinateSystem}, ip::ScalarInterpolation, dofs, qp::QuadraturePoint{ref_shape,T}, t) where {ref_shape,T}
+    x = @MVector zeros(T, 4)
+    @inbounds for i in 1:getnbasefunctions(ip)
+        val = Ferrite.shape_value(ip, qp.ξ, i)
+        x[1] += val * coeff.cs.u_transmural[dofs[i]]
+        x[2] += val * coeff.cs.u_apicobasal[dofs[i]]
+        x[3] += val * coeff.cs.u_rotational[dofs[i]]
+        x[4] += val * coeff.cs.u_transventricular[dofs[i]]
+    end
+    return BiVCoordinate(x[1], x[2], x[3], x[4])
+end
+
 """
     AnalyticalCoefficient(f::Function, cs::CoordinateSystemCoefficient)