GSOC-Task #2: Implement Faces and Edges recipies for Grid

[fverdugo]

In meeting on 2021-06-21, we agreed to introduce a new recipie, namely SurfaceWithEdges, surfacewithedges. Surface and wireframe plot will be a particular case.

I think that is better to have two separate recipies (this will be much more modular). I would call these recipies Faces and Edges for the plot types and faces and edges will be the functions automatically created by the @recipie macro.

With these recipies, we want to reproduce ALL plots in #16, in 1D, 2D, 3D, and different cell types, (e.g., quads and tris) in an unified API.

[fverdugo]

cc @ericneiva @amartinhuertas

[fverdugo]

This is the API we want to implement. The example is for 2D but the same would hold for 3D.

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
using CairoMakie
using GridapMakie
using FileIO

model = CartesianDiscreteModel((0.,1.5,0.,1.),(15,10)) |> simplexify
grid = get_grid(model)
celldata = rand(num_cells(grid))
nodaldata = rand(num_nodes(grid))

# Fig. 1: surface 
fig, = faces(grid)
save("fig_1.png",fig)

# Fig. 2: surface and edges 
fig, = faces(grid)
edges!(fig,grid)
save("fig_2.png",fig)

# Fig. 3: prescribe edge width 
fig, = faces(grid)
edges!(fig,grid,linewidth=2.5)
save("fig_3.png",fig)

# Fig. 4: set surface and edge color 
fig, = faces(grid,color=:green)
edges!(fig,grid,color=:red,linewidth=2.5)
save("fig_4.png",fig)

# Fig. 5: set surface color via a cell field
fig, = faces(grid,color=celldata,fieldstyle=:cells)
save("fig_5.png",fig)

# Fig. 6: set surface color via a nodal field
fig, = faces(grid,color=nodaldata,fieldstyle=:nodes)
save("fig_6.png",fig)

# Fig. 7: display colorbar 
fig,ax,tp = faces(grid,color=nodaldata)
Colorbar(fig[1,2],tp)
save("fig_7.png",fig)

# Fig. 8: change colormap
fig,ax,tp = faces(grid,color=nodaldata,colormap=:heat,colorrange=[0,1.])
Colorbar(fig[1,2],tp)
save("fig_8.png",fig)

# Fig. 9: edge color via nodal field
fig,ax,tp = edges(grid,color=nodaldata)
Colorbar(fig[1,2],tp)
save("fig_9.png",fig)

The generated figures should look like these ones:

Fig. 1: surface	Fig. 2: surface and edges	Fig. 3: prescribe edge width
Fig. 4: set surface and edge color	Fig. 5: set surface color via a cell field	Fig. 6: set surface color via a nodal field
Fig. 7: display colorbar	Fig. 8: change colormap	Fig. 9: edge color via nodal field

[fverdugo]

Some comments:

• The only new keyword argument I have considered is fieldstyle, with values :nodes (default) and :cells that specifies if the array provided in color has to be interpreted as nodal or cell data.

• The implementation of edges can rely on the implementation of faces for a grid of "edges". This will make the code very modular. That is edges(grid) will be implemented as an alias of faces(Grid(ReferenceFE{1},grid)) assuming that Grid(ReferenceFE{1},grid) returns the "edge" skeleton of grid (not sure if this method is already implemented in gridap, but Grid(ReferenceFE{1},model) should be working).

• At the end of the day, for 2D and 3D plots, we only need to draw triangles and segments embedded in 2D or 3D.

[fverdugo]

This function would extract the "edge" skeleton of a given grid:

function to_edge_grid(grid::UnstructuredGrid)
  topo = GridTopology(grid)
  labels = FaceLabeling(topo)
  model = DiscreteModel(grid,topo,labels)
  edge_grid = Grid(ReferenceFE{1},model)
end

thus using this function, edges(grid) will be an alias for faces(grid |> to_edge_grid)

[ericneiva]

I have edited the code snippet in #17 (comment) to fix some inconsistencies with the figure numbering. Now, I think the snippet and the figures below match.

[paurierap]

Comments on the current implementation:

• The type conversion from Gridap to GeometryBasics is no longer applied to all plot types PlotFunc, but just to Mesh. This is to allow our recipes to be defined for Grid types so that they may use certain functions (e.g. to_edge_grid).
• The dispatch in faces for an actual grid and the skeleton of a grid (the alias for edges(grid)) is given by the dimension of the cells Dc, with Dc=num_cell_dims(grid).
• The attributes of both faces and edges match those of Makie.mesh and Makie.wireframe. Additionally, the :colorrange key (needed to specify limits for a colorbar) is explicitly provided inside the plot! implementation and be given as an argument.
• The cycle key is set to nothing by default to avoid it overriding the default colors. This does not have negative consequences.
• Two new functions are defined in conversions.jl since they are used more than once: dimension_dispatchand get_nodes_and_ids.
• The implementation of faces for nodal and cell fields is different: for the former, the Makie.mesh function does all the work. However, to plot cell fields, we need to draw the grid simplex by simplex. Each simplex is associated with a color for a given colormap and is then plotted according to its coordinates. A future implementation will consider higher polytopes as a union of simplices (to be done soon).
• In edges, both Makie.linesegments seems to do the job (Makie.lines gives some continuity issues). Nonetheless, to achieve smoother corner caps of the edges, it is preferable to assign a color to every segment or face.

[paurierap]

Task 2 final API & Remarks

The implementation of both faces and edges is complete! Using these two recipes, we can draw 2D and 3D grids with the following API

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
using CairoMakie
using GridapMakie
using FileIO

model = CartesianDiscreteModel((0.,1.5,0.,1.),(15,10)) |> simplexify
grid = get_grid(model)
celldata = rand(num_cells(grid))
nodaldata = rand(num_nodes(grid))

# Fig. 1: surface 
fig = faces(grid)
save("fig_1.png", fig)

# Fig. 2: surface and edges 
fig, ax = faces(grid)
edges!(ax, grid)
save("fig_2.png", fig)

# Fig. 3: prescribe edge width 
fig, ax = faces(grid)
edges!(ax, grid, linewidth=2.5)
save("fig_3.png", fig)

# Fig. 4: set surface and edge color 
fig, ax = faces(grid, color=:lightseagreen)
edges!(ax, grid, color=:darkslategray)
save("fig_4.png", fig)

# Fig. 5: set surface color via a cell field
fig = faces(grid, color=celldata, fieldstyle=:cells)
save("fig_5.png", fig)

# Fig. 6: set surface color via a nodal field
fig = faces(grid, color=nodaldata, fieldstyle=:nodes)
save("fig_6.png", fig)

# Fig. 7: display colorbar 
fig, ax, plt = faces(grid, color=nodaldata, colorrange=(0,1))
Colorbar(fig[1,2], plt, ticks=0:0.25:1)
save("fig_7.png",fig)

# Fig. 8: change colormap
fig, ax, plt = faces(grid, color=nodaldata, colormap=:heat, colorrange=(0,1))
Colorbar(fig[1,2], plt, ticks=0:0.25:1)
save("fig_8.png",fig)

# Fig. 9: edge color via nodal field
fig, ax, plt = edges(grid, color=nodaldata, colorrange=(0,1))
Colorbar(fig[1,2], plt, ticks=0:0.25:1)
save("fig_9.png",fig)

The outcome of the figures will be:

Fig. 1: surface	Fig. 2: surface and edges	Fig. 3: prescribe edge width
Fig. 4: set surface and edge color	Fig. 5: set surface color via a cell field	Fig. 6: set surface color via a nodal field
Fig. 7: display colorbar	Fig. 8: change colormap	Fig. 9: edge color via nodal field

[amartinhuertas]

Great! Is the code in this branch ready to merge?

[paurierap]

Not yet. We still need to account for general polytopes in the cell field plot, which should be fixed as soon as we have the simplex-to-face mapping.

[ericneiva]

Excellent job, Pau! 👏 I am happy to see that you can do all the plots with the planned interface (with the only change of superposing the plots on the axis).

[fverdugo]

Wow both figures and API look really nice! Great Job!

[paurierap]

Pull requested! There is just the minor detail of cell fields with quads but should be soon fixed.

[fverdugo]

Hi @paurierap! Let's take a look at your PR tomorrow.

[fverdugo]

How to compute boundary grid and corresponding boundary face to cell map:

function to_boundary_grid_with_map(grid::Grid)
  Df = 2
  Dc = 3
  topo = GridTopology(grid)
  labels = FaceLabeling(topo)
  model = DiscreteModel(grid,topo,labels)
  face_grid = Grid(ReferenceFE{Df},model)
  face_to_mask = get_face_mask(labels,"boundary",Df)
  face_to_cells = get_faces(topo,Df,Dc)
  face_to_cell = lazy_map(first,face_to_cells)
  bface_to_face = findall(face_to_mask)
  bface_to_cell = lazy_map(Reindex(face_to_cell),bface_to_face)
  GridPortion(face_grid,bface_to_face), bface_to_cell
end

[fverdugo]

How to plot discontinuous fields (in particular cell-weise constant data)

using Gridap
using Gridap.Geometry
using Gridap.ReferenceFEs
using CairoMakie
using Makie.GeometryBasics
using FileIO

function to_point(x::VectorValue{D,T}) where {D,T}
  GeometryBasics.Point{D,T}(Tuple(x))
end

function to_dg_mesh(grid)
  ps = to_dg_points(grid)
  Dc = num_cell_dims(grid)
  Tc = Int32
  cns = Vector{Vector{Tc}}(undef,num_cells(grid))
  i = 1
  for cell in 1:num_cells(grid)
    cn = zeros(Tc,Dc+1)
    for lnode in 1:(Dc+1)
      cn[lnode] = i
      i += one(Tc)
    end
    cns[cell] = cn
  end
  fs = lazy_map(GeometryBasics.NgonFace{Dc+1,Tc},cns) |> collect
  GeometryBasics.connect(ps,fs) |> GeometryBasics.Mesh |> GeometryBasics.normal_mesh
end

function to_dg_points(grid)
  node_x = get_node_coordinates(grid)
  coords = to_dg_node_values(grid,node_x)
  ps = [ to_point(x) for x in coords]
  ps
end

function to_dg_node_values(grid,node_value)
  cell_nodes = get_cell_node_ids(grid)
  i = 0
  cache = array_cache(cell_nodes)
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache,cell_nodes,cell)
    for node in nodes
      i += 1
    end
  end
  T = eltype(node_value)
  values = zeros(T,i)
  i = 0
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache,cell_nodes,cell)
    for node in nodes
      i += 1
      values[i] = node_value[node]
    end
  end
  values
end

function to_dg_cell_values(grid,cell_value)
  cell_nodes = get_cell_node_ids(grid)
  i = 0
  cache = array_cache(cell_nodes)
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache,cell_nodes,cell)
    for node in nodes
      i += 1
    end
  end
  T = eltype(cell_value)
  values = zeros(T,i)
  i = 0
  for cell in 1:length(cell_nodes)
    nodes = getindex!(cache,cell_nodes,cell)
    for node in nodes
      i += 1
      values[i] = cell_value[cell]
    end
  end
  values
end

n=5
domain = (0,1,0,1)
cells = (n,n)
model = CartesianDiscreteModel(domain,cells) |> simplexify

grid2 = Grid(ReferenceFE{2},model)
grid1 = Grid(ReferenceFE{1},model)
grid0 = Grid(ReferenceFE{0},model)

data2 = 2*rand(num_faces(model,2))
data1 = 2*rand(num_faces(model,1))
data0 = 2*rand(num_faces(model,0))

mesh2 = to_dg_mesh(grid2)
values2 = to_dg_node_values(grid2,data0)
fig,ax,sc = mesh(mesh2,color=values2)
save("f_2d_2_0.png",fig)

values2 = to_dg_cell_values(grid2,data2)
fig,ax,sc = mesh(mesh2,color=values2)
save("f_2d_2_2.png",fig)

ps1 = to_dg_points(grid1)
values1 = to_dg_node_values(grid1,data0)
fig,ax,sc = linesegments(ps1,color=values1)
save("f_2d_1_0.png",fig)

ps1 = to_dg_points(grid1)
values1 = to_dg_cell_values(grid1,data1)
fig,ax,sc = linesegments(ps1,color=values1)
save("f_2d_1_1.png",fig)

ps0 = to_dg_points(grid0)
values0 = to_dg_node_values(grid0,data0)
fig,ax,sc = scatter(ps0,color=values0)
save("f_2d_0_0.png",fig)

fig,ax,sc = mesh(mesh2,color=rand(length(coordinates(mesh2))))
linesegments!(ps1,color=:yellow)
scatter!(ps0,color=:cyan)
save("f_2d_012.png",fig)

# now in 3d

domain = (0,1,0,1,0,1)
cells = (n,n,n)
model = CartesianDiscreteModel(domain,cells) |> simplexify

grid2 = Grid(ReferenceFE{2},model)
grid1 = Grid(ReferenceFE{1},model)
grid0 = Grid(ReferenceFE{0},model)

data2 = 2*rand(num_faces(model,2))
data1 = 2*rand(num_faces(model,1))
data0 = 2*rand(num_faces(model,0))

mesh2 = to_dg_mesh(grid2)
values2 = to_dg_node_values(grid2,data0)
fig,ax,sc = mesh(mesh2,color=values2)
save("f_3d_2_0.png",fig)

values2 = to_dg_cell_values(grid2,data2)
fig,ax,sc = mesh(mesh2,color=values2)
save("f_3d_2_2.png",fig)

ps1 = to_dg_points(grid1)
values1 = to_dg_node_values(grid1,data0)
fig,ax,sc = linesegments(ps1,color=values1)
save("f_3d_1_0.png",fig)

ps1 = to_dg_points(grid1)
values1 = to_dg_cell_values(grid1,data1)
fig,ax,sc = linesegments(ps1,color=values1)
save("f_3d_1_1.png",fig)

ps0 = to_dg_points(grid0)
values0 = to_dg_node_values(grid0,data0)
fig,ax,sc = scatter(ps0,color=values0)
save("f_3d_0_0.png",fig)

fig,ax,sc = mesh(mesh2,color=rand(length(coordinates(mesh2))))
linesegments!(ps1,color=:yellow,linewidth=4)
scatter!(ps0,color=:cyan,markersize=30)
save("f_3d_012.png",fig)

Edit: Added normal_mesh at the end of to_dg_mesh to account for 3D grids.

[fverdugo]

• use to_dg_mesh to plot both nodal and cell data (this would simplify the code)
• Rename dimension_dispatch -> to_plot_grid (use type parameter to dispatch among 2D and 3D)
• add vertices recipe. It should be an alias for faces(grid |> to_vertex_grid), and define to_vertex_grid which should be very similar to to_edge_grid.
• Add function simplexify_with_map in Gridap. So that it returns also the mapping from simplex cells to cube cells. In any case, we agreed to handle n-cubes in a later step.

[fverdugo]

• Rename faces -> cells.
• Get rid of fieldstyle attribute.

[paurierap]

Comments on the current implementation v2:

• There is no type conversion anymore. Instead, we transform our grid via to_dg_mesh to a suitable GeometryBasics object to be plotted by, respectively in 0D, 1D and 2D (3D grids are reduced to 2D faces), Makie.scatter, Makie.linesegments and Makie.mesh.
• The dispatch in cells for a grid and the skeleton of a grid (either 0D or 1D, i.e. the alias edges(grid) or vertices(grid)) is given by the dimension of the cells Dc, with Dc=num_cell_dims(grid).
• All the attributes of cells, edges and vertices match those of Makie.Mesh, Makie.Lines and Makie.Scatter. For example, linewidth modifies the width of the lines drawn by edges or markersize determines the diameter in pixels of the dots in scatter.
• The cycle attribute is set to nothing by default to avoid it overriding the default colors. The same with shading, which appears to be set to true by default. Finally, colorrange=nothing so that it can be easily modified (either automatically or by the user).
• The functions to_cell_grid, to_edge_grid and to_vertex_grid allow to extract the 2D, 1D and 0D skeletons of a grid.
• Lastly, the color features are handled by color_handling, which accounts for the type of plot that we are looking for depending on the length of color::Vector (or even if color::Symbol or color=RGBAf0(1,0.2,1,1), for example). Also, we got rid of the fieldstyle attribute, that we thought was confusing.
• Important: Currently we cannot display colorbars for colormaps drawn by edges unless the colorrange attribute is explicitly defined. This is because in the edges recipe, we ought to set colorrange to nothing (the same would happen for vertices). As for the third task we only consider single color edges, we can postpone the solution of this issue.

[fverdugo]

Fixed by PR #24