Hey there,
I am using a dolfinx script to simulate a port hamiltonian system. Based on example scripts from Mr. Dokken’s github page, I am trying to visualize my results of a beam (over time) using pyvista, where I copied from the example the function
def harry_plotter(t_input, u_input):
"""
Create a figure of the concentration u_input warped visualized in 3D at timet step t.
"""
plotter = pyvista.Plotter()
# Create grid defined by the function space for visualization of the function
topology, cells, geometry = dlf.plot.create_vtk_mesh(u_input.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
var_name = f"u({t_input})"
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_input)] = u_input.x.array.reshape(geometry.shape[0], len(u_input))
function_grid[var_name] = values
function_grid.set_active_vectors(var_name)
# Warp mesh by deformation
warped = function_grid.warp_by_vector(var_name, factor=1)
# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped)
plotter.show_axes()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
pyvista.start_xvfb()
figure_as_array = plotter.screenshot(f"diffusion_{t_input:.2f}.png")
# Clear plotter for next plot
plotter.remove_actor(actor)
I just renamed some of the variables due to their names already existing in my script.
Since I defined my elements as
P_elem = ufl.VectorElement("CG", gebiet.ufl_cell(), 1)
R_elem = ufl.TensorElement("DG", gebiet.ufl_cell(), 1)
P_wish = ufl.VectorElement("CG", gebiet.ufl_cell(), 1)
R_wish = ufl.TensorElement("DG", gebiet.ufl_cell(), 1)
P_wish_2 = ufl.VectorElement("CG", gebiet.ufl_cell(), 1)
mel = ufl.MixedElement([P_elem, R_elem, P_wish, R_wish, P_wish_2])
U = dlf.fem.FunctionSpace(gebiet, mel)
u = dlf.fem.Function(U)
(p, F, v, P, pos) = ufl.split(u)
(dv, dP, dp, dF, dpos) = ufl.TestFunctions(U)
u_n = dlf.fem.Function(U)
(p_n, F_n, v_n, P_n, pos_n) = ufl.split(u_n),
where gebiet is the domain and pos is the position of the elements in space, I call the plotting function with harry_plotter(t_temp, u.sub(0)). The full script, if necessary for helping me, is below this post. Some of the variables are in German to avoid accidentally using a name already in use.
This unfortunately does not work, as I am getting the error:
/usr/bin/python3 hyperelasticity2.py
Zeit: 0.0
Traceback (most recent call last):
File "/home/marco/Schreibtisch/Code
/Kernproblem/hyperelasticity2.py", line 178, in <module>
harry_plotter(t_temp, u.sub(0))
File "/home/marco/Schreibtisch/Code
/Kernproblem/hyperelasticity2.py", line 121, in harry_plotter
topology, cells, geometry = dlf.plot.create_vtk_mesh(u_input.function_space)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/plot.py", line 105, in _
dofmap_ = dofmap.list.array.reshape(dofmap.list.num_nodes, num_dofs_per_cell)
ValueError: cannot reshape array of size 51840 into shape (4320,4)
I am really confused. Inputting u in total doesn’t work, as expected, but all the u.sub(...) do not work either, where I would have expected to work at least the vector elements, especially the position, u.sub(4).
Do you have any idea how to make this plotting work?
Thank you so much in advance.
Best regards
Marco
Full script:
## import stuff
import dolfinx as dlf
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt
import ufl
import scipy.io
from petsc4py import PETSc
import pyvista
## define constants
L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lmbda = beta
g = gamma
rho_0=1
mu=1
lamb=1
nu=np.array([1,1,1])
delta_t = 1/16
t_gesamt = 1/4 # +delta_t/2 # s
N_t = round(t_gesamt / delta_t)
k_w = 3.0
rho_w = 2.0
## define mesh
gebiet = dlf.mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],[20,6,6])
# gebiet = dlf.mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],[20,6,6], cell_type=dlf.mesh.CellType.hexahedron)
Normalenvektor = ufl.FacetNormal(gebiet)
## define function spaces
P_elem = ufl.VectorElement("CG", gebiet.ufl_cell(), 1)
R_elem = ufl.TensorElement("DG", gebiet.ufl_cell(), 1)
P_wish = ufl.VectorElement("CG", gebiet.ufl_cell(), 1)
R_wish = ufl.TensorElement("DG", gebiet.ufl_cell(), 1)
P_wish_2 = ufl.VectorElement("CG", gebiet.ufl_cell(), 1)
mel = ufl.MixedElement([P_elem, R_elem, P_wish, R_wish, P_wish_2])
U = dlf.fem.FunctionSpace(gebiet, mel)
u = dlf.fem.Function(U)
(p, F, v, P, pos) = ufl.split(u)
(dv, dP, dp, dF, dpos) = ufl.TestFunctions(U)
u_n = dlf.fem.Function(U)
(p_n, F_n, v_n, P_n, pos_n) = ufl.split(u_n)
V_tau = dlf.fem.FunctionSpace(gebiet,P_elem)
tau_0 = dlf.fem.Function(V_tau) # Spannungsvektor
## define initial conditions
def einervektor(x):
values = np.zeros((gebiet.geometry.dim, x.shape[1]), dtype=np.float64)
values[0,:] = 1
values[1,:] = 1
values[2,:] = 1
return values
tau_0.interpolate(einervektor)
def initial_cond(x):
values = np.zeros((gebiet.geometry.dim*gebiet.geometry.dim, x.shape[1]), dtype=np.float64)
values[0,:] = 1
values[4,:] = 1
values[8,:] = 1
return values
u.sub(1).interpolate(initial_cond)
u_n.sub(1).interpolate(initial_cond)
## define functions
Psi_fun = lambda F: (mu / 2) * (ufl.tr(ufl.dot(F.T,F)) - 3) - mu * ufl.ln(ufl.sqrt(ufl.classes.Determinant(ufl.dot(F.T,F)))) + (lmbda / 2) * (ufl.ln(ufl.sqrt(ufl.classes.Determinant(ufl.dot(F.T,F)))))**2
P_fun = lambda F: ufl.dot(F,mu*(ufl.classes.Identity(F.ufl_shape[0])-ufl.inv(ufl.dot(F.T,F)))+lamb*ufl.ln(ufl.sqrt(ufl.classes.Determinant(ufl.dot(F.T,F))))*ufl.inv(ufl.dot(F.T,F)))
# psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# P = ufl.diff(psi, F)
## define BCs
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], L)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], W)),
(5, lambda x: np.isclose(x[2], 0)),
(6, lambda x: np.isclose(x[2], W))] # 2 weitere Ränder für 3. Dimension, gleicher Index für selbes Randelement oder Elemente mit logischem Oder verknüpfen.
facet_indices, facet_markers = [], []
fdim = gebiet.topology.dim - 1
for (marker, locator) in boundaries:
facets = dlf.mesh.locate_entities(gebiet, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dlf.mesh.meshtags(gebiet, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=gebiet, subdomain_data=facet_tag)
ds_Dirichlet = ds(0)
ds_Neumann = ds(1)
dx=ufl.Measure('dx', domain=gebiet, metadata={'quadrature_degree': 4})
# fdim = domain.topology.dim -1
# left_facets = mesh.locate_entities_boundary(domain, fdim, left)
# right_facets = mesh.locate_entities_boundary(domain, fdim, right)
def harry_plotter(t_input, u_input):
"""
Create a figure of the concentration u_input warped visualized in 3D at timet step t.
"""
plotter = pyvista.Plotter()
# Create grid defined by the function space for visualization of the function
topology, cells, geometry = dlf.plot.create_vtk_mesh(u_input.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
var_name = f"u({t_input})"
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u_input)] = u_input.x.array.reshape(geometry.shape[0], len(u_input))
function_grid[var_name] = values
function_grid.set_active_vectors(var_name)
# Warp mesh by deformation
warped = function_grid.warp_by_vector(var_name, factor=1)
# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped)
plotter.show_axes()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
pyvista.start_xvfb()
figure_as_array = plotter.screenshot(f"diffusion_{t_input:.2f}.png")
# Clear plotter for next plot
plotter.remove_actor(actor)
## loop over time steps
time = np.zeros(N_t + 1)
Hamiltonian = np.zeros(N_t + 1)
for count, t_temp in enumerate(np.linspace(0, t_gesamt, N_t + 1)):
print('Zeit: ',t_temp)
## update load vector
# skipped
## define residuum
F_FEM = ufl.dot(dv, (p-p_n) / delta_t) * dx - ufl.inner(ufl.grad(dv), P) * dx - ufl.dot(dv, ufl.dot(P,Normalenvektor)) * ds_Neumann + ufl.dot(dv,tau_0) * ds_Neumann + \
ufl.inner(dP, (F-F_n) / delta_t) * dx + ufl.inner(dP, ufl.grad(v)) * dx - ufl.dot(ufl.dot(dP,Normalenvektor),v) * ds_Dirichlet + \
ufl.dot(dp,v-p/rho_0) * ufl.dx + \
ufl.inner(dF,P-P_fun(F)) * dx + \
ufl.dot(dpos, (pos - pos_n)/ delta_t) * dx - ufl.dot(dpos, v) * ufl.dx # integration der Position, Zeitintegration noch auf implizite Mittelpunktsregel umstellen! # nu+
## solve problem
problem = dlf.fem.petsc.NonlinearProblem(F_FEM, u)
solver = dlf.nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
n, converged = solver.solve(u) # Wo ist der Lösungsoutput? u wird mit dem Ergebnis überschrieben.
## Update old field
u_n.x.array[:] = u.x.array[:]
(p_n, F_n, v_n, P_n, pos_n) = ufl.split(u_n)
## update Hamiltonian
test_konstant = dlf.fem.Constant(gebiet, 1/2)
time[count] = t_temp
# https://fenicsproject.discourse.group/t/incompatible-function-arguments-for-assemble-scalar/8194 -> https://github.com/FEniCS/dolfinx/blob/7b3e22e95f62f827a4e9d3eb5124f53c5e2148f3/python/test/unit/fem/test_assemble_submesh.py#L55
Hamiltonian[count] = dlf.fem.assemble_scalar(dlf.fem.form(1 / 2 * ufl.dot(p_n, p_n) / rho_0 * dx + Psi_fun(F_n) * dx))
## Plot beam
harry_plotter(t_temp, u.sub(0))
## plot hamiltonian over time
plt.figure()
plt.title('Hamiltonian', fontsize=10)
plt.plot(time, Hamiltonian)
plt.xlabel("$t$")
plt.ylabel("$H$")
plt.grid(True)
#plt.savefig('results/Hamiltonian.pdf')
plt.show()