Fixes to domain decomposition, non-scalar variable updates

Project.toml

@@ -1,7 +1,7 @@
 name = "Jutul"
 uuid = "2b460a1a-8a2b-45b2-b125-b5c536396eb9"
 authors = ["Olav Møyner <[email protected]>"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -41,6 +41,7 @@ Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
+Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
 KaHyPar = "2a6221f6-aa48-11e9-3542-2d9e0ef01880"
 LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e"
@@ -54,6 +55,7 @@ PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 JutulAMGCLWrapExt = "AMGCLWrap"
 JutulGLMakieExt = "GLMakie"
 JutulGraphMakieExt = ["GraphMakie", "NetworkLayout", "LayeredLayouts", "Makie"]
+JutulGmshExt = "Gmsh"
 JutulHYPREExt = "HYPRE"
 JutulKaHyParExt = "KaHyPar"
 JutulMakieExt = "Makie"
@@ -101,6 +103,7 @@ julia = "1.7"
 CPUSummary = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
+Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
 KaHyPar = "2a6221f6-aa48-11e9-3542-2d9e0ef01880"
 LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e"

ext/JutulGmshExt/JutulGmshExt.jl

@@ -0,0 +1,14 @@
+module JutulGmshExt
+    using Jutul
+    import Jutul: UnstructuredMesh, IndirectionMap, check_equal_perm
+    using Gmsh: Gmsh, gmsh
+
+    using StaticArrays
+    using OrderedCollections
+
+    const QUAD_T = SVector{4, Int}
+    const TRI_T = SVector{3, Int}
+
+    include("utils.jl")
+    include("interface.jl")
+end

ext/JutulGmshExt/interface.jl

@@ -0,0 +1,66 @@
+
+function Jutul.mesh_from_gmsh(pth; verbose = false, kwarg...)
+    Gmsh.initialize()
+    ext = pth |> splitext |> last
+    gmsh.open(pth)
+    if lowercase(ext) == ".geo"
+        gmsh.model.mesh.generate()
+    end
+
+    dim = gmsh.model.getDimension()
+    dim == 3 || error("Only 3D models are supported")
+
+    gmsh.model.mesh.removeDuplicateNodes()
+    node_tags, pts, = gmsh.model.mesh.getNodes()
+    node_remap = Dict{UInt64, Int}()
+    for (i, tag) in enumerate(node_tags)
+        tag::UInt64
+        node_remap[tag] = i
+    end
+    remaps = (
+        nodes = node_remap,
+        faces = Dict{UInt64, Int}(),
+        cells = Dict{UInt64, Int}()
+    )
+    pts = reshape(pts, Int(dim), :)
+    pts_s = collect(vec(reinterpret(SVector{3, Float64}, pts)))
+
+    @assert size(pts, 2) == length(node_tags)
+    faces_to_nodes = parse_faces(remaps, verbose = verbose)
+    face_lookup = generate_face_lookup(faces_to_nodes)
+
+    cells_to_faces = parse_cells(remaps, faces_to_nodes, face_lookup, verbose = verbose)
+    # We are done with Gmsh.jl at this point, finalize it.
+    Gmsh.finalize()
+
+    neighbors = build_neighbors(cells_to_faces, faces_to_nodes, face_lookup)
+
+    # Make both of these in case we have rogue faces that are not connected to any cell.
+    bnd_faces = Int[]
+    int_faces = Int[]
+    n_bad = 0
+    for i in axes(neighbors, 2)
+        l, r = neighbors[:, i]
+        l_bnd = l == 0
+        r_bnd = r == 0
+
+        if l_bnd || r_bnd 
+            if l_bnd && r_bnd
+                n_bad += 1
+            else
+                push!(bnd_faces, i)
+            end
+        else
+            push!(int_faces, i)
+        end
+    end
+    bnd_neighbors, bnd_faces_to_nodes, bnd_cells_to_faces = split_boundary(neighbors, faces_to_nodes, cells_to_faces, bnd_faces, boundary = true)
+    int_neighbors, int_faces_to_nodes, int_cells_to_faces = split_boundary(neighbors, faces_to_nodes, cells_to_faces, int_faces, boundary = false)
+
+    c2f = IndirectionMap(int_cells_to_faces)
+    c2b = IndirectionMap(bnd_cells_to_faces)
+    f2n = IndirectionMap(int_faces_to_nodes)
+    b2n = IndirectionMap(bnd_faces_to_nodes)
+    print_message("Mesh parsed successfully:\n    $(length(c2f)) cells\n    $(length(f2n)) internal faces\n    $(length(b2n)) boundary faces\n    $(length(pts_s)) nodes", verbose)
+    return UnstructuredMesh(c2f, c2b, f2n, b2n, pts_s, int_neighbors, bnd_neighbors; kwarg...)
+end

ext/JutulGmshExt/utils.jl

@@ -0,0 +1,225 @@
+
+function add_next!(faces, remap, tags, numpts, offset)
+    vals = Int[]
+    for j in 1:numpts
+        push!(vals, remap[tags[offset + j]])
+    end
+    push!(faces, vals)
+end
+
+function parse_faces(remaps; verbose = false)
+    node_remap = remaps.nodes
+    face_remap = remaps.faces
+    faces = Vector{Int}[]
+    for (dim, tag) in gmsh.model.getEntities()
+        if dim != 2
+            continue
+        end
+        type = gmsh.model.getType(dim, tag)
+        name = gmsh.model.getEntityName(dim, tag)
+        elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
+        type = gmsh.model.getType(dim, tag)
+        ename = gmsh.model.getEntityName(dim, tag)
+        for (etypes, etags, enodetags) in zip(elemTypes, elemTags, elemNodeTags)
+            name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(etypes)
+            if name == "Quadrilateral 4"
+                numpts = 4
+            elseif name == "Triangle 3"
+                numpts = 3
+            else
+                error("Unsupported element type $name for faces.")
+            end
+            @assert length(enodetags) == numpts*length(etags)
+            print_message("Faces: Processing $(length(etags)) tags of type $name", verbose)
+            for (i, etag) in enumerate(etags)
+                offset = (i-1)*numpts
+                add_next!(faces, node_remap, enodetags, numpts, offset)
+                face_remap[etag] = length(faces)
+            end
+            print_message("Added $(length(etags)) faces of type $name with $(length(unique(enodetags))) unique nodes", verbose)
+        end
+    end
+    return faces
+end
+
+function get_cell_decomposition(name)
+    if name == "Hexahedron 8"
+        tris = Tuple{}()
+        quads = (
+            QUAD_T(0, 4, 7, 3),
+            QUAD_T(1, 2, 6, 5),
+            QUAD_T(0, 1, 5, 4),
+            QUAD_T(2, 3, 7, 6),
+            QUAD_T(0, 3, 2, 1),
+            QUAD_T(4, 5, 6, 7)
+        )
+        numpts = 8
+    elseif name == "Tetrahedron 4"
+        tris = (
+            TRI_T(0, 1, 3),
+            TRI_T(0, 2, 1),
+            TRI_T(0, 3, 2),
+            TRI_T(1, 2, 3)
+        )
+        quads = Tuple{}()
+        numpts = 4
+    elseif name == "Pyramid 5"
+        # TODO: Not really tested.
+        tris = (
+            TRI_T(0, 1, 4),
+            TRI_T(0, 4, 3),
+            TRI_T(3, 4, 2),
+            TRI_T(1, 2, 4)
+        )
+        quads = (QUAD_T(0, 3, 2, 1),)
+        numpts = 4
+    elseif name == "Prism 6"
+        # TODO: Not really tested.
+        tris = (
+            TRI_T(0, 2, 1),
+            TRI_T(3, 4, 5)
+        )
+        quads = (
+            QUAD_T(0, 1, 4, 3),
+            QUAD_T(0, 3, 5, 2),
+            QUAD_T(1, 2, 5, 4)
+        )
+        numpts = 6
+    else
+        error("Unsupported element type $name for cells.")
+    end
+    return (tris, quads, numpts)
+end
+
+function print_message(msg, verbose)
+    if verbose
+        println(msg)
+    end
+end
+
+function parse_cells(remaps, faces, face_lookup; verbose = false)
+    node_remap = remaps.nodes
+    face_remap = remaps.faces
+    cell_remap = remaps.cells
+    cells = Vector{Tuple{Int, Int}}[]
+    for (dim, tag) in gmsh.model.getEntities()
+        if dim != 3
+            continue
+        end
+        type = gmsh.model.getType(dim, tag)
+        name = gmsh.model.getEntityName(dim, tag)
+        # Get the mesh elements for the entity (dim, tag):
+        elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag)
+        # * Type and name of the entity:
+        type = gmsh.model.getType(dim, tag)
+        ename = gmsh.model.getEntityName(dim, tag)
+        for (etypes, etags, enodetags) in zip(elemTypes, elemTags, elemNodeTags)
+            name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(etypes)
+            tris, quads, numpts = get_cell_decomposition(name)
+            print_message("Cells: Processing $(length(etags)) tags of type $name", verbose)
+            @assert length(enodetags) == numpts*length(etags)
+            nadded = 0
+            for (i, etag) in enumerate(etags)
+                offset = (i-1)*numpts
+                pt_range = (offset+1):(offset+numpts)
+                @assert length(pt_range) == numpts
+                pts = map(i -> node_remap[enodetags[i]], pt_range)
+                cell = Tuple{Int, Int}[]
+                cellno = length(cells)
+                for face_t in (tris, quads)
+                    for (fno, face) in enumerate(face_t)
+                        face_pts = map(i -> pts[i+1], face)
+                        face_pts_sorted = sort(face_pts)
+                        faceno = get(face_lookup, face_pts_sorted, 0)
+                        if faceno == 0
+                            nadded += 1
+                            push!(faces, face_pts)
+                            faceno = length(faces)
+                            face_lookup[face_pts_sorted] = faceno
+                            sgn = 1
+                        else
+                            sgn = check_equal_perm(face_pts, faces[faceno]) ? 1 : 2
+                        end
+                        push!(cell, (faceno, sgn))
+                    end
+                end
+                push!(cells, cell)
+            end
+            print_message("Added $(length(etags)) new cells of type $name and $nadded new faces.", verbose)
+        end
+    end
+    return cells
+end
+
+function build_neighbors(cells, faces, face_lookup)
+    neighbors = zeros(Int, 2, length(faces))
+    for (cellno, cell_to_faces) in enumerate(cells)
+        for (face_index, lr) in cell_to_faces
+            face_pts = faces[face_index]
+            face_pts_sorted = sort(face_pts)
+            faceno = face_lookup[face_pts_sorted]
+            face_ref = faces[faceno]
+            oldn = neighbors[lr, faceno]
+            oldn == 0 || error("Cannot overwrite face neighbor for cell $cellno - was already defined as $oldn for index $lr: $(neighbors[:, faceno])")
+            neighbors[lr, faceno] = cellno
+        end
+    end
+    return neighbors
+end
+
+function generate_face_lookup(faces)
+    face_lookup = Dict{Union{QUAD_T, TRI_T}, Int}()
+
+    for (i, face) in enumerate(faces)
+        n = length(face)
+        if n == 3
+            ft = sort(TRI_T(face[1], face[2], face[3]))
+        elseif n == 4
+            ft = sort(QUAD_T(face[1], face[2], face[3], face[4]))
+        else
+            error("Unsupported face type with $n nodes, only 3 (for tri) and 4 (for quad) are known.")
+        end
+        @assert issorted(ft)
+        face_lookup[ft] = i
+    end
+    return face_lookup
+end
+
+function split_boundary(neighbors, faces_to_nodes, cells_to_faces, active_ix::Vector{Int}; boundary::Bool)
+    remap = OrderedDict{Int, Int}()
+    for (i, ix) in enumerate(active_ix)
+        remap[ix] = i
+    end
+    # is_active = [false for _ in eachindex(faces_to_nodes)]
+    # is_active[active_ix] .= true
+    # Make renumbering here.
+    if boundary
+        new_neighbors = Int[]
+        for ix in active_ix
+            l, r = neighbors[:, ix]
+            @assert l == 0 || r == 0
+            push!(new_neighbors, max(l, r))
+        end
+    else
+        new_neighbors = Tuple{Int, Int}[]
+        for ix in active_ix
+            l, r = neighbors[:, ix]
+            @assert l != 0 && r != 0
+            push!(new_neighbors, (l, r))
+        end
+    end
+    new_faces_to_nodes = map(copy, faces_to_nodes[active_ix])
+    # Handle cells -> current type of faces
+    new_cells_to_faces = Vector{Int}[]
+    for cell_to_faces in cells_to_faces
+        new_cell = Int[]
+        for (face, sgn) in cell_to_faces
+            if haskey(remap, face)
+                push!(new_cell, remap[face])
+            end
+        end
+        push!(new_cells_to_faces, new_cell)
+    end
+
+    return (new_neighbors, new_faces_to_nodes, new_cells_to_faces)
+end

ext/JutulMakieExt/mesh_plots.jl

@@ -27,14 +27,28 @@ function Jutul.plot_mesh_impl!(ax, m;
     has_face_filter = !isnothing(faces)
     has_bface_filter = !isnothing(boundaryfaces)
     if has_cell_filter || has_face_filter || has_bface_filter
+        keep_cells = Dict{Int, Bool}()
+        keep_faces = Dict{Int, Bool}()
+        keep_bf = Dict{Int, Bool}()
+
         if eltype(cells) == Bool
             @assert length(cells) == number_of_cells(m)
             cells = findall(cells)
         end
+        if has_cell_filter
+            for c in cells
+                keep_cells[c] = true
+            end
+        end
         if eltype(faces) == Bool
             @assert length(faces) == number_of_faces(m)
             faces = findall(faces)
         end
+        if has_face_filter
+            for f in faces
+                keep_faces[f] = true
+            end
+        end
         if eltype(boundaryfaces) == Bool
             @assert length(boundaryfaces) == number_of_boundary_faces(m)
             boundaryfaces = findall(boundaryfaces)
@@ -43,6 +57,9 @@ function Jutul.plot_mesh_impl!(ax, m;
             nf = number_of_faces(m)
             boundaryfaces = deepcopy(boundaryfaces)
             boundaryfaces .+= nf
+            for f in boundaryfaces
+                keep_bf[f] = true
+            end
         end
         ntri = size(tri, 1)
         keep = fill(false, ntri)
@@ -53,13 +70,13 @@ function Jutul.plot_mesh_impl!(ax, m;
             tri_tmp = tri[i, 1]
             keep_this = true
             if has_cell_filter
-                keep_this = keep_this && cell_ix[tri_tmp] in cells
+                keep_this = keep_this && haskey(keep_cells, cell_ix[tri_tmp])
             end
             if has_face_filter
-                keep_this = keep_this && face_ix[tri_tmp] in faces
+                keep_this = keep_this && haskey(keep_faces, face_ix[tri_tmp])
             end
             if has_bface_filter
-                keep_this = keep_this && face_ix[tri_tmp] in boundaryfaces
+                keep_this = keep_this && haskey(keep_bf, face_ix[tri_tmp])
             end
             keep[i] = keep_this
         end