Visualization of geometries deformed by displacement fields

[fverdugo]

Simple code sniped of how to deform your visualization grid with a generic displacement field.

using Gridap
using Gridap.Geometry
using Gridap.ReferenceFEs
using GridapMakie
using GLMakie
using Gridap.Visualization

domain = (0,1,0,1,0,1)
cell_nums = 10 .* (1,1,1)
model = CartesianDiscreteModel(domain,cell_nums)
Ω = Triangulation(model)

α  = π*0.25
r = TensorValue(cos(α),-sin(α),0., sin(α),cos(α),0., 0.,0.,1.)

fun(x)=r⋅x
uh = CellField(fun,Ω)

vd = visualization_data(Ω,"",cellfields=["u"=>uh])
grid = vd[1].grid
vals = vd[1].nodaldata["u"]
display(get_node_coordinates(grid))
display(vals)

grid2 = deepcopy(grid)
nodes2 = get_node_coordinates(grid2)
nodes2 .= nodes2 .+ vals

fig = plot(grid2)
display(fig)

We want to implement a nice high level API for this. Something linke:

Ω2 = warp(Ω,uh)
plot(Ω2)

cc @Kevin-Mattheus-Moerman

[Kevin-Mattheus-Moerman]

@fverdugo thanks for following up with this.

FYI I tried to build something to show a deformed model. If I build a model like this:

# Model
domain = (0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)

# Setup integration
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

Then size(Ω.node_coords[:]) produces (64,), i.e. there are 64 nodes for a mesh defining the corners of the hexahedra.

If I do the following however then size(nodes2) is (216,)

vd = visualization_data(Ω,"",cellfields=["u"=>uh])
grid = vd[1].grid
vals = vd[1].nodaldata["u"]
grid2 = deepcopy(grid)
nodes2 = get_node_coordinates(grid2)
nodes2 .= nodes2 .+ vals

Is there a way to obtain displacements and coordinates in the form (64,) ? Alternatively, how can I get the cell descriptions for the quadratic hex elements?

[Kevin-Mattheus-Moerman]

Got it. If I get element descriptions like this, then the indices conform to the (216,) system:

  vd=visualization_data(Ω,"");
  grid = vd[1].grid
  E = get_cell_node_ids(grid) #Elements
  V = get_node_coordinates(grid) #Nodes

[fverdugo]

how can I get the cell descriptions for the quadratic hex elements?

Try this:

vd=visualization_data(Ω,"",order=2);

[Kevin-Mattheus-Moerman]

@fverdugo at first I thought it was the order but actually Ω.node_coords gives the shared nodal coordinate set (64 nodes in for a 3x3x3 element = 4x4x4 shared node mesh), while for the approach with visualization_data the nodes are "unshared" (leading to a set of (3x3x3)x8=216 nodes).

[fverdugo]

Yes, and the general one is the "unshared" since it also allows one to visualise discontinuous functions as needed in advanced discretization methods like DG.

[Kevin-Mattheus-Moerman]

Thanks @fverdugo. What would I use to get visualization data for the shared node case? Or do we currently default to unshared?