How to properly get coordinates with a ::Grid structure

New to Julia
1

Dear all,

I’m a newbie at Julia. Probably I’m doing something wrong and/or inefficiently: I was trying to replicate one of the examples to generate my own linear elasticity solver. Starting with a very simple example (following incompressible elasticity example ).

function create_cook_grid(nx, ny)
    corners = [Vec{2}((0.0,   0.0)),
               Vec{2}((48.0, 44.0)),
               Vec{2}((48.0, 60.0)),
               Vec{2}((0.0,  44.0))]
    grid = generate_grid(Triangle, (nx, ny), corners);
    # facesets for boundary conditions
    addfaceset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);
    addfaceset!(grid, "traction", x -> norm(x[1]) ≈ 48.0);
    return grid
end;

and copying assembly functions, I can found a solution using P1 and P2 finite elements. In VTK looks great. However, when I try to get nodal coordinates and its connectivity:

grid.nodes; # has 99-element Array{Node{2,Float64},1}
grid.cells;   # has 160-element Array{Cell{2,3,3},1}

while the length of the 2d vector solution, called u, is 5202. The latter coincides with the info provided by Paraview: num of cells = 5000, num of points = 2601.

Trying to collect the nodal coordinates (certainly this looks sub-optimal):

nodes = grid.nodes
elems = grid.cells;
    
nnodes = length(nodes); # << 99-element Array{Node{2,Float64},1}
nelems = length(elems); # << 160-element Array{Cell{2,3,3},1}
    
coords = Array{Float64,2}(undef,nnodes,2) 
cells  = Array{Int,2}(undef,nelems,3)
    
for i in 1:nnodes
   coords[i,:] = [nodes[i].x[1],nodes[i].x[2]];
end
    
for i in 1:nelems
   cells[i,:] = [elems[i].nodes[1],elems[i].nodes[2],elems[i].nodes[3]]
end

and then

using PyPlot
PyPlot.plot_trisurf(coords[:,1],coords[:,2],triangles=cells)

shows a nice mesh, but coarser than that shown in Paraview.

I’m completely lost. How can I get the real nodal coordinates and connectivity?
Another question: Had someone import the VTK library in Julia? I’m thinking of having something similar to Mayavi in Julia.

Thanks.
J

2

I’m sorry, I found myself the answer searching along the sources:

np = getnnodes(grid) # 2601
nc = getncells(grid) # 5000
coord1 = getcoordinates(grid,1) # Array{Tensor{1,2,Float64,2},1} with coordinates of the node 1
coords = getnodes(grid) # 2601-element Array{Node{2,Float64},1} : with all the set of nodes
cells =  getcells(grid) # 5000-element Array{Cell{2,3,3},1}: with all the cells (elements)

So it is consistent now.

Anyway, I was wondering why it is beneficial to have a Finite Element dataset composed by Arrays of Node or Cells structures, which in turn are Tensor structures, instead of merely Array of Array structures. Is it faster in this way? Easier to parallelize? Maybe I’m thinking of old-fashioned C++ codes.

Finally, what about plotting meshes and solutions 2d and 3d in Jupyter? Using Makie or Plotly is a bit cumbersome since we need to pass Arrays instead of Nodes and Tensors. It would be great to have some ideas or recommendations about that.

Thanks.