add trivial.jl solving the trivial equation f = u

[jw3126]

Aimed at the documentation part of https://github.com/gridap/Gridap.jl/issues/345

[fverdugo]

Dear @jw3126, thanks for your PR.

We believe that it would be very nice to have a tutorial on how to print a PDE solution in pure Julia (i.e., without using Paraview). But we would include it once we can visualize at least in 2D.

If you like, you can help us to explore how we can print 1D, 2D, and possibly 3D solutions in pure Julia, perhaps using https://github.com/JuliaPlots/Makie.jl

It would be very helpful if you can contribute to this.

On the other hand, the main focus of the tutorial would be visualization. The FE machinery is already presented e.g. here https://gridap.github.io/Tutorials/dev/pages/t001_poisson/ So, we don't need to consider a trivial equation f=u (which is perhaps too trivial). If you want a trivial equation, we better consider the Poisson equation.

[jw3126]

Yeah, it is tricky to choose a nice set of tutorials. To me, it is often useful in any context to see the most trivial example. But I also see your point about redundancy. What I dislike about the poisson tutorial is, that it is not copy pastable. Because you need the mesh file.

[fverdugo]

BTW, in https://gridap.github.io/Tutorials/stable/pages/t002_validation/ we consider a Poisson eq on a Cartesian grid.

In any case, the part we like from your PR is to try to visualize things in pure Julia. The particular PDE is not so important.

[jw3126]

I will look into makie. I think linear Lagrange elements should be not so difficult.

[fverdugo]

I will look into makie. I think linear Lagrange elements should be not so difficult.

If you find how to visualize on top of Lagrangian elements, this is all what we need. In Gridap, we convert any other element types to Lagrangian before visualizing.

Start with this example:

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
domain = (0,1,0,1)
cells = (10,10)
grid = CartesianGrid(domain,cells) # For quadrilaterals
# grid = simplexify(CartesianGrid(domain,cells)) # For triangles
x = get_node_coordinates(trian)
f(x) = sin(4*x[1])*cos(3*x[2])
fx = apply(f,x) 

The first step is to visualize fx on top of grid with makie.

Once you are able to solve this problem, I can tell you how to couple this with the Gridap visualization machinery.

[jw3126]

I tried the following. The missing piece is the connectivity. E.g. which node indices belong to a common cell. How to extract that from a grid?

using Makie

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
domain = (0,1,0,1)
cells = (1,1)
model = CartesianDiscreteModel(domain,cells)
trian = Triangulation(model)

pts = get_node_coordinates(model)
f(pt) = pt[2]
makie_coords = [x[i][j] for i in 1:length(x), j in 1:2]
makie_conn = [
    1 2 3;
    3 4 2;
]
# color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
color = vec((map(f, pts)))
@assert size(makie_coords) == (length(color), 2)
@assert makie_coords isa Matrix{Float64}
@assert color isa Vector{Float64}
scene = mesh(makie_coords, makie_conn, color = color, shading = false)
wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 3)
display(scene)

[jw3126]

A quick and drity way is to use meshscatter.

using Makie

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
domain = (0,1,0,1)
cells = (100,100)
grid = CartesianDiscreteModel(domain,cells) # For quadrilaterals
# grid = simplexify(CartesianGrid(domain,cells)) # For triangles
trian = Triangulation(grid)
pts = get_node_coordinates(trian)
f(x) = sin(4*x[1])*cos(3*x[2])
x = collect(vec(apply(pt -> pt[1], pts)))
y = collect(vec(apply(pt -> pt[2], pts)))
z = collect(vec(apply(f,pts) ))
scene = meshscatter(x,y,z, color=z, markersize=1e-2)
display(scene)

But it looks horrible for coarse grids.

[fverdugo]

The missing piece is the connectivity. E.g. which node indices belong to a common cell. How to extract that from a grid?

get_cell_nodes(trian)

[jw3126]

using Makie

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry

function to_makie_matrix(T, itr)
    x1 = first(itr)
    out = Matrix{T}(undef, length(itr), length(x1))
    for i in 1:length(itr)
        for j in 1:length(x1)
            out[i, j] = itr[i][j]
        end
    end
    out
end

domain = (0,2pi,0,pi)
cells = (10,20)
model = CartesianDiscreteModel(domain,cells)
model = simplexify(model)
trian = Triangulation(model)

pts = get_node_coordinates(model)
f(pt) = sin(pt[1]) * cos(pt[2])
makie_coords = to_makie_matrix(Float64, pts)
makie_conn = to_makie_matrix(Int, get_cell_nodes(trian))
color = vec(map(f, pts))
scene = mesh(makie_coords, makie_conn, color = color, shading = false)
wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 3)
display(scene)

makie_coords = hcat(makie_coords, color)

[fverdugo]

Great! does it work for squares? and what about for 3D?

[jw3126]

It does neither work for squares nor 3d. What would you expect for 3d? The meshscatter would work in 3d.

[jw3126]

Still I think this would be a useful tutorial, short term. Mid term it would be nice to have a GridapMakie package.

[fverdugo]

Yes, perhaps this is enough at this moment for the Makie side.

However, some (very minor) work is still to be done in the Gridap side in order to be able to plot any kind of FE function, not only Lagrangian ones.

Perhaps you can also help with this. It shouldn't be very difficult. I can explain how to proceed.

[fverdugo]

Mid term it would be nice to have a GridapMakie package.

Yes, this is the final goal.

[jw3126]

However, some (very minor) work is still to be done in the Gridap side in order to be able to plot any kind of FE function, not only Lagrangian ones.

Even for linear lagrangian, I am not yet fully sure how to do it. How do I extract the function values at the node coordinates?