Meshes and mesh utilities

Mesh types

Jutul has two main internal mesh types: Cartesian meshes and unstructured meshes. The unstructured format is a general polyhedral mesh format, and a Cartesian mesh can easily be converted to an unstructured mesh. Coarsened meshes can be created by a fine scale mesh and a partition vector.

Jutul.CartesianMesh — Type
CartesianMesh(dims, [Δ, [origin]])

Create a Cartesian mesh with dimensions specified by the Tuple dims.

Arguments

  • dims::Tuple: Number of grid cells in each direction. For example, (nx, ny) will give a 2D grids with nx cells in the x-direction.
  • Δ::Tuple=Tuple(ones(length(dims))): Equal length to dims. First option: A

Tuple of scalars where each entry is the length of each cell in that direction. For example, specifying (Δx, Δy) for a uniform grid with each grid cell having area ofΔx*Δy. Second option: ATuple` of vectors where each entry contains the cell sizes in the direction.

  • origin=zeros(length(dims)): The origin of the first corner in the grid.

Examples

Generate a uniform 3D mesh that discretizes a domain of 2 by 3 by 5 units with 3 by 5 by 2 cells:

julia> CartesianMesh((3, 5, 2), (2.0, 3.0, 5.0))
CartesianMesh (3D) with 3x5x2=30 cells

Generate a non-uniform 2D mesh:

julia> CartesianMesh((2, 3), ([1.0, 2.0], [0.1, 3.0, 2.5]))
CartesianMesh (2D) with 2x3x1=6 cells
source
Jutul.UnstructuredMesh — Type
UnstructuredMesh(g::CartesianMesh)

Convert CartesianMesh instance to unstructured grid. Note that the mesh must be 2D and 3D for a 1-to-1 conversion. 1D meshes are implicitly converted to 2D.

source
Jutul.CoarseMesh — Type
CoarseMesh(G::JutulMesh, p)

Construct a coarse mesh from a given JutulMesh that can be converted to an UnstructuredMesh instance. The second argument p should be a partition Vector with one entry per cell in the original grid that assigns that cell to a coarse block. Should be one-indexed and the numbering should be sequential and contain at least one fine cell for each coarse index. This is tested by the function.

source
Jutul.MRSTWrapMesh — Type
MRSTWrapMesh(G, N = nothing)

Mesh that adapts an exported MRST mesh to the Jutul interface. G is assumed to be read directly from file using MAT.matread. The raw exported grid can be found under the data field.

source

Plotting functions

Plotting requires that a Makie backend is loaded (typically GLMakie or CairoMakie). The documentation uses CairoMakie to work on machines without OpenGL enabled, but if you want fast and interactive plots, GLMakie should be preferred.

Non-mutating

Jutul.plot_mesh — Function
plot_mesh(mesh)
plot_mesh(mesh;
    cells = nothing,
    faces = nothing,
    boundaryfaces = nothing,
    outer = false,
    color = :lightblue,
)

Plot a mesh with uniform colors. Optionally, indices cells, faces or boundaryfaces can be passed to limit the plotting to a specific selection of entities.

source
Jutul.plot_cell_data — Function
plot_cell_data(mesh::JutulMesh, data::Vector; kwarg...)
plot_cell_data(mesh, data;
    cells = nothing,
    faces = nothing,
    boundaryfaces = nothing
)

Plot cell-wise values (as a vector) on the mesh. Optionally, indices cells, faces or boundaryfaces can be passed to limit the plotting to a specific selection of entities.

source
Jutul.plot_interactive — Function
plot_interactive(mesh, vector_of_dicts; kwarg...)

Launch an interactive plot of a mesh with the given vector_of_dicts (or just a dict). Each dict can have cell data either as vectors (one value per cell) or matrices (one column per cell).

source

Mutating

Jutul.plot_mesh! — Function
plot_mesh!(ax, mesh)

Mutating version of plot_mesh that plots into an existing Makie Axis instance.

source
Jutul.plot_cell_data! — Function
plot_cell_data!(ax, mesh, data; kwarg...)

Mutating version of plot_cell_data that plots into an existing Makie Axis

source
Jutul.plot_mesh_edges! — Function
plot_mesh_edges!(ax, mesh; kwarg...)

Plot the edges of all cells on the exterior of a mesh into existing Makie Axis ax.

source

Example: Cartesian meshes

For example, we can make a small 2D mesh with given physical dimensions and convert it:

using Jutul
nx = 10
ny = 5
g2d_cart = CartesianMesh((nx, ny), (100.0, 50.0))
CartesianMesh (2D) with 10x5x1=50 cells
g2d = UnstructuredMesh(g2d_cart)
UnstructuredMesh with 50 cells, 85 faces and 30 boundary faces

We can then plot it, colorizing each cell by its enumeration:

using CairoMakie
fig, ax, plt = plot_cell_data(g2d, 1:number_of_cells(g2d))
plot_mesh_edges!(ax, g2d)
fig
Example block output

If we want to drill down a bit further, we can make a plot:

We can make a 3D mesh in the same manner:

nz = 3
g3d = UnstructuredMesh(CartesianMesh((nx, ny, nz), (100.0, 50.0, 30.0)))
UnstructuredMesh with 150 cells, 355 faces and 190 boundary faces

And plot it the same way:

using CairoMakie
nc = number_of_cells(g3d)
fig, ax, plt = plot_cell_data(g3d, 1:nc)
plot_mesh_edges!(ax, g3d)
fig
Example block output

We can also plot only a subset of cells:

using CairoMakie
fig, ax, plt = plot_cell_data(g3d, 1:nc, cells = 1:2:nc)
fig
Example block output

Mesh API functions

Queries

Jutul.number_of_cells — Function
number_of_cells(D::Union{DataDomain, DiscretizedDomain})

Get the number of cells in a DataDomain or DiscretizedDomain.

source
number_of_cells(g)::Integer

Get the number of cells in a mesh.

source
Jutul.number_of_faces — Function
number_of_faces(D::Union{DataDomain, DiscretizedDomain})

Get the number of faces in a DataDomain or DiscretizedDomain.

source
number_of_faces(g)::Integer

Get the number of faces in a mesh.

source

Manipulation

Jutul.extrude_mesh — Function
extrude_mesh(m2d::UnstructuredMesh, nlayers)
extrude_mesh(m2d::UnstructuredMesh, [1, 2, 5, 10])

Extrude a 2D mesh into a 3D mesh by adding layers of cells in the z-direction. The number of layers can be specified as an integer or as an array of depths. The depths must be in increasing order.

source
Jutul.extract_submesh — Function
extract_submesh(g::UnstructuredMesh, cells)

Extract a subgrid for a given mesh and a iterable of cells to keep.

source

Example: Mesh manipulation

We can quickly build new meshes by applying transformations to an already existing mesh. Let us create a Cartesian mesh and extract the cells that lie within a circle:

using Jutul, CairoMakie
g = UnstructuredMesh(CartesianMesh((10, 10), (1.0, 1.0)))
geo = tpfv_geometry(g)
keep = Int[]
for c in 1:number_of_cells(g)
    x, y = geo.cell_centroids[:, c]
    if (x - 0.5)^2 + (y - 0.5)^2 < 0.25
        push!(keep, c)
    end
end
subg = extract_submesh(g, keep)
fig, ax, plt = plot_mesh(subg)
plot_mesh_edges!(ax, g, color = :red)
fig
Example block output

We can turn this into a 3D mesh by extruding it, and then tweak the nodes:

g3d = Jutul.extrude_mesh(subg, 20)
for node in eachindex(subg.node_points)
    g3d.node_points[node] += 0.01*rand(3)
end
fig, ax, plt = plot_mesh(g3d)
fig
Example block output

Geometry

Jutul.TwoPointFiniteVolumeGeometry — Type
TwoPointFiniteVolumeGeometry(neighbors, areas, volumes, normals, cell_centers, face_centers)

Store two-point geometry information for a given list of neighbors specified as a 2 by n matrix where n is the number of faces such that face i connectes cells N[1, i] and N[2, i].

The two-point finite-volume geometry contains the minimal set of geometry information required to compute standard finite-volume discretizations.

source
Jutul.find_enclosing_cells — Function
find_enclosing_cells(G, traj; geometry = tpfv_geometry(G), n = 25)

Find the cell indices of cells in the mesh G that are intersected by a given trajectory traj. traj can be either a matrix with equal number of columns as dimensions in G (i.e. three columns for 3D) or a Vector of SVector instances with the same length.

The optional argument geometry is used to define the centroids and normals used in the calculations. You can precompute this if you need to perform many searches. The keyword argument n can be used to set the number of discretizations in each segment.

use_boundary is by default set to false. If set to true, the boundary faces of cells are treated more rigorously when picking exactly what cells are cut by a trajectory, but this requires that the boundary normals are oriented outwards, which is currently not the case for all meshes from downstream packages.

limit_box speeds up the search by limiting the search to the minimal bounding box that contains both the trajectory and the mesh. This can be turned off by passing false. There should be no difference in the cells tagged by changing this option.

source

Example: Cell intersection

using CairoMakie, Jutul
# 3D mesh
G = CartesianMesh((4, 4, 5), (100.0, 100.0, 100.0))
trajectory = [
    50.0 25.0 1;
    55 35.0 25;
    65.0 40.0 50.0;
    70.0 70.0 90.0
]

cells = Jutul.find_enclosing_cells(G, trajectory)

# Optional plotting, requires Makie:
fig, ax, plt = Jutul.plot_mesh_edges(G)
plot_mesh!(ax, G, cells = cells, alpha = 0.5, transparency = true)
lines!(ax, trajectory, linewidth = 10)
fig
Example block output

2D version:

using CairoMakie, Jutul
# 2D mesh
G = CartesianMesh((50, 50), (1.0, 2.0))
trajectory = [
    0.1 0.1;
    0.2 0.4;
    0.3 1.2
]
fig, ax, plt = Jutul.plot_mesh_edges(G)
cells = Jutul.find_enclosing_cells(G, trajectory)
# Plotting, needs Makie
fig, ax, plt = Jutul.plot_mesh_edges(G)
plot_mesh!(ax, G, cells = cells, alpha = 0.5, transparency = true)
lines!(ax, trajectory[:, 1], trajectory[:, 2], linewidth = 3)
fig
Example block output

Mesh generation

Gmsh support

Jutul.mesh_from_gmsh — Function
G = mesh_from_gmsh(pth)
G = mesh_from_gmsh()
G = mesh_from_gmsh(pth; verbose = true)

Parse a Gmsh file and return a Jutul UnstructuredMesh (in 3D only). Requires the Gmsh.jl package to be loaded. If no path is provided in pth it is assumed that you are managing the Gmsh state manually and it will use the current selected mesh inside Gmsh. Please note that Gmsh is GPL licensed unless you have obtained another type of license from the authors.

source

Radial mesh

using Jutul, CairoMakie
import Jutul.RadialMeshes: radial_mesh
import Jutul: plot_mesh_edges
nangle = 10
radii = [0.2, 0.5, 1.0]
m = radial_mesh(nangle, radii; centerpoint = true)
plot_mesh_edges(m)
(Scene (1600px, 900px):
  0 Plots
  1 Child Scene:
    └ Scene (1600px, 900px), Makie.Axis3(), MakieCore.LineSegments{Tuple{Base.ReinterpretArray{GeometryBasics.Point{2, Float64}, 1, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}, Vector{Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}}, false}}})

Radial meshes can still be indexed as Cartesian meshes:

IJ = map(i -> cell_ijk(m, i), 1:number_of_cells(m))
30-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 1, 1)
 (2, 1, 1)
 (3, 1, 1)
 (4, 1, 1)
 (5, 1, 1)
 (6, 1, 1)
 (7, 1, 1)
 (8, 1, 1)
 (9, 1, 1)
 (10, 1, 1)
 ⋮
 (2, 3, 1)
 (3, 3, 1)
 (4, 3, 1)
 (5, 3, 1)
 (6, 3, 1)
 (7, 3, 1)
 (8, 3, 1)
 (9, 3, 1)
 (10, 3, 1)
fig, ax, plt = plot_cell_data(m, map(first, IJ))
fig
Example block output
fig, ax, plt = plot_cell_data(m, map(ijk -> ijk[2], IJ))
fig
Example block output

We can also plot the faces by using the nodes together with standard Makie plotting calls. We regenerate the mesh and make it contain a single cell in the middle before plotting it:

m = radial_mesh(nangle, radii; centerpoint = false)
ncells = number_of_cells(m)
fig, ax, plt = plot_cell_data(m, 1:ncells, alpha = 0.25)
scatter!(ax, m.node_points)
for face in 1:number_of_faces(m)
    n1, n2 = m.faces.faces_to_nodes[face]
    pt1 = m.node_points[n1]
    pt2 = m.node_points[n2]
    lines!(ax, [pt1, pt2], color = :red)
end

for bface in 1:number_of_boundary_faces(m)
    n1, n2 = m.boundary_faces.faces_to_nodes[bface]
    pt1 = m.node_points[n1]
    pt2 = m.node_points[n2]
    lines!(ax, [pt1, pt2], color = :blue)
end
fig
Example block output

We can also zoom in on a single cell and plot the oriented normals:

using LinearAlgebra
geo = tpfv_geometry(m)
cellno = 1
fig, ax, plt = plot_mesh(m, cells = cellno)
for face in m.faces.cells_to_faces[cellno]
    n1, n2 = m.faces.faces_to_nodes[face]
    @info "Interior edge $n1 to $n2"
    pt1 = m.node_points[n1]
    pt2 = m.node_points[n2]
    lines!(ax, [pt1, pt2], color = :red)
    midpt = (pt1 + pt2) / 2
    if m.faces.neighbors[face][1] == cellno
        sgn = 1
    else
        sgn = -1
    end
    lines!(ax, [midpt, midpt + sgn*norm(pt2 - pt1, 2)*geo.normals[:, face]], color = :orange)
end
for bface in m.boundary_faces.cells_to_faces[cellno]
    n1, n2 = m.boundary_faces.faces_to_nodes[bface]
    @info "Exterior edge $n1 to $n2"
    pt1 = m.node_points[n1]
    pt2 = m.node_points[n2]
    lines!(ax, [pt1, pt2], color = :blue)
    midpt = (pt1 + pt2) / 2
    lines!(ax, [midpt, midpt + norm(pt2 - pt1, 2)*geo.boundary_normals[:, bface]], color = :green)
end
scatter!(ax, m.node_points)
fig
Example block output

Spiral meshes

Jutul.RadialMeshes.spiral_mesh — Function
mesh = spiral_mesh(10, 3, spacing = [0.0, 0.5, 1.0])

Spiral mesh generator. Generates a spiral mesh in 2D with an optional "spacing" subdiscretization between each segment.

source
using Jutul, CairoMakie
import Jutul.RadialMeshes: spiral_mesh
n_angular_sections = 10
nrotations = 4
spacing = [0.0, 0.5, 1.0]
rmesh = spiral_mesh(n_angular_sections, nrotations, spacing = spacing)
num_cells = number_of_cells(rmesh)

fig, ax, plt = plot_cell_data(rmesh, 1:number_of_cells(rmesh))
fig
Example block output
import Jutul.RadialMeshes: spiral_mesh_tags
tags = spiral_mesh_tags(rmesh, spacing)
fig = Figure(size = (400, 1800))
for (figno, pp) in enumerate(pairs(tags))
    k, val = pp
    ax = Axis(fig[figno, 1], title = "Spiral tag $k")
    plot_cell_data!(ax, rmesh, val)
end
fig
Example block output