|
| 1 | + |
| 2 | +function add_next!(faces, remap, tags, numpts, offset) |
| 3 | + vals = Int[] |
| 4 | + for j in 1:numpts |
| 5 | + push!(vals, remap[tags[offset + j]]) |
| 6 | + end |
| 7 | + push!(faces, vals) |
| 8 | +end |
| 9 | + |
| 10 | +function parse_faces(remaps; verbose = false) |
| 11 | + node_remap = remaps.nodes |
| 12 | + face_remap = remaps.faces |
| 13 | + faces = Vector{Int}[] |
| 14 | + for (dim, tag) in gmsh.model.getEntities() |
| 15 | + if dim != 2 |
| 16 | + continue |
| 17 | + end |
| 18 | + type = gmsh.model.getType(dim, tag) |
| 19 | + name = gmsh.model.getEntityName(dim, tag) |
| 20 | + elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag) |
| 21 | + type = gmsh.model.getType(dim, tag) |
| 22 | + ename = gmsh.model.getEntityName(dim, tag) |
| 23 | + for (etypes, etags, enodetags) in zip(elemTypes, elemTags, elemNodeTags) |
| 24 | + name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(etypes) |
| 25 | + if name == "Quadrilateral 4" |
| 26 | + numpts = 4 |
| 27 | + elseif name == "Triangle 3" |
| 28 | + numpts = 3 |
| 29 | + else |
| 30 | + error("Unsupported element type $name for faces.") |
| 31 | + end |
| 32 | + @assert length(enodetags) == numpts*length(etags) |
| 33 | + print_message("Faces: Processing $(length(etags)) tags of type $name", verbose) |
| 34 | + for (i, etag) in enumerate(etags) |
| 35 | + offset = (i-1)*numpts |
| 36 | + add_next!(faces, node_remap, enodetags, numpts, offset) |
| 37 | + face_remap[etag] = length(faces) |
| 38 | + end |
| 39 | + print_message("Added $(length(etags)) faces of type $name with $(length(unique(enodetags))) unique nodes", verbose) |
| 40 | + end |
| 41 | + end |
| 42 | + return faces |
| 43 | +end |
| 44 | + |
| 45 | +function get_cell_decomposition(name) |
| 46 | + if name == "Hexahedron 8" |
| 47 | + tris = Tuple{}() |
| 48 | + quads = ( |
| 49 | + QUAD_T(0, 4, 7, 3), |
| 50 | + QUAD_T(1, 2, 6, 5), |
| 51 | + QUAD_T(0, 1, 5, 4), |
| 52 | + QUAD_T(2, 3, 7, 6), |
| 53 | + QUAD_T(0, 3, 2, 1), |
| 54 | + QUAD_T(4, 5, 6, 7) |
| 55 | + ) |
| 56 | + numpts = 8 |
| 57 | + elseif name == "Tetrahedron 4" |
| 58 | + tris = ( |
| 59 | + TRI_T(0, 1, 3), |
| 60 | + TRI_T(0, 2, 1), |
| 61 | + TRI_T(0, 3, 2), |
| 62 | + TRI_T(1, 2, 3) |
| 63 | + ) |
| 64 | + quads = Tuple{}() |
| 65 | + numpts = 4 |
| 66 | + elseif name == "Pyramid 5" |
| 67 | + # TODO: Not really tested. |
| 68 | + tris = ( |
| 69 | + TRI_T(0, 1, 4), |
| 70 | + TRI_T(0, 4, 3), |
| 71 | + TRI_T(3, 4, 2), |
| 72 | + TRI_T(1, 2, 4) |
| 73 | + ) |
| 74 | + quads = (QUAD_T(0, 3, 2, 1),) |
| 75 | + numpts = 4 |
| 76 | + elseif name == "Prism 6" |
| 77 | + # TODO: Not really tested. |
| 78 | + tris = ( |
| 79 | + TRI_T(0, 2, 1), |
| 80 | + TRI_T(3, 4, 5) |
| 81 | + ) |
| 82 | + quads = ( |
| 83 | + QUAD_T(0, 1, 4, 3), |
| 84 | + QUAD_T(0, 3, 5, 2), |
| 85 | + QUAD_T(1, 2, 5, 4) |
| 86 | + ) |
| 87 | + numpts = 6 |
| 88 | + else |
| 89 | + error("Unsupported element type $name for cells.") |
| 90 | + end |
| 91 | + return (tris, quads, numpts) |
| 92 | +end |
| 93 | + |
| 94 | +function print_message(msg, verbose) |
| 95 | + if verbose |
| 96 | + println(msg) |
| 97 | + end |
| 98 | +end |
| 99 | + |
| 100 | +function parse_cells(remaps, faces, face_lookup; verbose = false) |
| 101 | + node_remap = remaps.nodes |
| 102 | + face_remap = remaps.faces |
| 103 | + cell_remap = remaps.cells |
| 104 | + cells = Vector{Tuple{Int, Int}}[] |
| 105 | + for (dim, tag) in gmsh.model.getEntities() |
| 106 | + if dim != 3 |
| 107 | + continue |
| 108 | + end |
| 109 | + type = gmsh.model.getType(dim, tag) |
| 110 | + name = gmsh.model.getEntityName(dim, tag) |
| 111 | + # Get the mesh elements for the entity (dim, tag): |
| 112 | + elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag) |
| 113 | + # * Type and name of the entity: |
| 114 | + type = gmsh.model.getType(dim, tag) |
| 115 | + ename = gmsh.model.getEntityName(dim, tag) |
| 116 | + for (etypes, etags, enodetags) in zip(elemTypes, elemTags, elemNodeTags) |
| 117 | + name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(etypes) |
| 118 | + tris, quads, numpts = get_cell_decomposition(name) |
| 119 | + print_message("Cells: Processing $(length(etags)) tags of type $name", verbose) |
| 120 | + @assert length(enodetags) == numpts*length(etags) |
| 121 | + nadded = 0 |
| 122 | + for (i, etag) in enumerate(etags) |
| 123 | + offset = (i-1)*numpts |
| 124 | + pt_range = (offset+1):(offset+numpts) |
| 125 | + @assert length(pt_range) == numpts |
| 126 | + pts = map(i -> node_remap[enodetags[i]], pt_range) |
| 127 | + cell = Tuple{Int, Int}[] |
| 128 | + cellno = length(cells) |
| 129 | + for face_t in (tris, quads) |
| 130 | + for (fno, face) in enumerate(face_t) |
| 131 | + face_pts = map(i -> pts[i+1], face) |
| 132 | + face_pts_sorted = sort(face_pts) |
| 133 | + faceno = get(face_lookup, face_pts_sorted, 0) |
| 134 | + if faceno == 0 |
| 135 | + nadded += 1 |
| 136 | + push!(faces, face_pts) |
| 137 | + faceno = length(faces) |
| 138 | + face_lookup[face_pts_sorted] = faceno |
| 139 | + sgn = 1 |
| 140 | + else |
| 141 | + sgn = check_equal_perm(face_pts, faces[faceno]) ? 1 : 2 |
| 142 | + end |
| 143 | + push!(cell, (faceno, sgn)) |
| 144 | + end |
| 145 | + end |
| 146 | + push!(cells, cell) |
| 147 | + end |
| 148 | + print_message("Added $(length(etags)) new cells of type $name and $nadded new faces.", verbose) |
| 149 | + end |
| 150 | + end |
| 151 | + return cells |
| 152 | +end |
| 153 | + |
| 154 | +function build_neighbors(cells, faces, face_lookup) |
| 155 | + neighbors = zeros(Int, 2, length(faces)) |
| 156 | + for (cellno, cell_to_faces) in enumerate(cells) |
| 157 | + for (face_index, lr) in cell_to_faces |
| 158 | + face_pts = faces[face_index] |
| 159 | + face_pts_sorted = sort(face_pts) |
| 160 | + faceno = face_lookup[face_pts_sorted] |
| 161 | + face_ref = faces[faceno] |
| 162 | + oldn = neighbors[lr, faceno] |
| 163 | + oldn == 0 || error("Cannot overwrite face neighbor for cell $cellno - was already defined as $oldn for index $lr: $(neighbors[:, faceno])") |
| 164 | + neighbors[lr, faceno] = cellno |
| 165 | + end |
| 166 | + end |
| 167 | + return neighbors |
| 168 | +end |
| 169 | + |
| 170 | +function generate_face_lookup(faces) |
| 171 | + face_lookup = Dict{Union{QUAD_T, TRI_T}, Int}() |
| 172 | + |
| 173 | + for (i, face) in enumerate(faces) |
| 174 | + n = length(face) |
| 175 | + if n == 3 |
| 176 | + ft = sort(TRI_T(face[1], face[2], face[3])) |
| 177 | + elseif n == 4 |
| 178 | + ft = sort(QUAD_T(face[1], face[2], face[3], face[4])) |
| 179 | + else |
| 180 | + error("Unsupported face type with $n nodes, only 3 (for tri) and 4 (for quad) are known.") |
| 181 | + end |
| 182 | + @assert issorted(ft) |
| 183 | + face_lookup[ft] = i |
| 184 | + end |
| 185 | + return face_lookup |
| 186 | +end |
| 187 | + |
| 188 | +function split_boundary(neighbors, faces_to_nodes, cells_to_faces, active_ix::Vector{Int}; boundary::Bool) |
| 189 | + remap = OrderedDict{Int, Int}() |
| 190 | + for (i, ix) in enumerate(active_ix) |
| 191 | + remap[ix] = i |
| 192 | + end |
| 193 | + # is_active = [false for _ in eachindex(faces_to_nodes)] |
| 194 | + # is_active[active_ix] .= true |
| 195 | + # Make renumbering here. |
| 196 | + if boundary |
| 197 | + new_neighbors = Int[] |
| 198 | + for ix in active_ix |
| 199 | + l, r = neighbors[:, ix] |
| 200 | + @assert l == 0 || r == 0 |
| 201 | + push!(new_neighbors, max(l, r)) |
| 202 | + end |
| 203 | + else |
| 204 | + new_neighbors = Tuple{Int, Int}[] |
| 205 | + for ix in active_ix |
| 206 | + l, r = neighbors[:, ix] |
| 207 | + @assert l != 0 && r != 0 |
| 208 | + push!(new_neighbors, (l, r)) |
| 209 | + end |
| 210 | + end |
| 211 | + new_faces_to_nodes = map(copy, faces_to_nodes[active_ix]) |
| 212 | + # Handle cells -> current type of faces |
| 213 | + new_cells_to_faces = Vector{Int}[] |
| 214 | + for cell_to_faces in cells_to_faces |
| 215 | + new_cell = Int[] |
| 216 | + for (face, sgn) in cell_to_faces |
| 217 | + if haskey(remap, face) |
| 218 | + push!(new_cell, remap[face]) |
| 219 | + end |
| 220 | + end |
| 221 | + push!(new_cells_to_faces, new_cell) |
| 222 | + end |
| 223 | + |
| 224 | + return (new_neighbors, new_faces_to_nodes, new_cells_to_faces) |
| 225 | +end |
0 commit comments