Problem when exporting displacemen data regular grid from SLEPc modal analysis using boxfield


57
views
0
8 weeks ago by
Daniel  
I want to export structured data from a 3D modal analysis in linear elasticity. For setting up the modal analysis I got inspired by https://comet-fenics.readthedocs.io/en/latest/ The computed mode shape I would like to export on a regular grid. As suggested in the FEniCS tutorial I used the boxfield module which can be found in the github repository for the tutorial. Additionally, I export the data as a pvd file for ParaView. The source code looks like
import fenics as fe
import numpy as np
import boxfield
import os

# Numerical parameters
n_eigen = 10
i_eigen = 0

# Geometry parameters
l = 2500.0e-6
b = 1270.0e-6
h = 20.0e-6

# Mesh parameters
nx = 100
ny = 50
nz = 2

# Material parameters
E = 130e9
nu = 0.28
rho = 2330.0

# Lame parameters
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))

# Definitions of stress and strain
def epsilon(u):
    return fe.sym(fe.grad(u))

def sigma(u):
    n_dim = u.geometric_dimension()
    return 2 * mu * epsilon(u) + lmbda * fe.tr(epsilon(u)) * fe.Identity(n_dim)

# Create mesh
mesh = fe.BoxMesh(fe.Point(0, 0, 0), fe.Point(l, b, h), nx, ny, nz)

# Create test and trial functions
V = fe.VectorFunctionSpace(mesh, 'Lagrange', degree=1)
u = fe.TrialFunction(V)
v = fe.TestFunction(V)

# Define one clamped end as boundary condition
def boundary(x, on_boundary):
    return fe.near(x[0], 0)

bc = fe.DirichletBC(V, fe.Constant((0.0, 0.0, 0.0)), boundary)

# Assemble stiffness matrix
k_form = fe.inner(sigma(u), epsilon(v)) * fe.dx
l_form = fe.Constant(1) * u[0] * fe.dx
K = fe.PETScMatrix()
b_tmp = fe.PETScVector()
fe.assemble_system(k_form, l_form, bc, A_tensor=K, b_tensor=b_tmp)

# Assemble mass matrix
m_form = rho * fe.dot(u, v) * fe.dx
M = fe.PETScMatrix()
fe.assemble(m_form, tensor=M)

# Setup SLEPc eigensolver
eigensolver = fe.SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectrum'] = 'smallest real'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.0

# Compute eigenvalues
eigensolver.solve(n_eigen)

# Pick on eigenvalue
r, c, rx, cx = eigensolver.get_eigenpair(i_eigen)

# Determine eigenmode in fenics
u_ex = fe.Function(V)
u_ex.vector()[:] = rx

# Export eigenmode as structured np using boxfield
u_box = boxfield.FEniCSBoxField(u_ex, (nx, ny, nz))
x = u_box.grid.coor[boxfield.X]
y = u_box.grid.coor[boxfield.Y]
z = u_box.grid.coor[boxfield.Z]
uu = u_box.values[0, :, :, 1]
vv = u_box.values[1, :, :, 1]
ww = u_box.values[2, :, :, 1]
np.savez('mode.npz', x=x, y=y, z=z, u=uu, v=vv, w=ww)

# Expoert eigenmode pvd file
f_pvd = fe.File('mode.pvd')
f_pvd << u_ex​

In ParaView the mode shape looks as expected for a cantilevered structure.

The white structure is the original mesh and the bend structure is the computed mode shape which looks like expected for a simple cantilevered body.

However, when I plot the numpy data I obtain these colormaps showing the different vector components of the deflection.

Note the non-continuous deflection pattern not present in the ParaView image. This becomes even more clear when looking at a cross section of the displacement data along the x direction. The data does not even fulfill the boundary condition anymore and has a periodic pattern with the deflection jumping to values close to zero.

The code for loading and plotting the data is

import numpy as np
import matplotlib.pyplot as plt

# Load data
data = np.load('mode.npz')
x = data['x']
y = data['y']
z = data['z']
u = data['u']
v = data['v']
w = data['w']
data.close()

# Normalize data
u_max = np.max(np.abs(u))
v_max = np.max(np.abs(v))
w_max = np.max(np.abs(w))
uvw_max = np.max([u_max, v_max, w_max])

u = u / uvw_max
v = v / uvw_max
w = w / uvw_max

# Displacement colormaps
plt.figure(1, figsize=(11, 5))
plt.clf()

plt.subplot(1, 3, 1)
plt.title('Displacement in x direction')
plt.imshow(u)
plt.colorbar()
plt.xlabel('y')
plt.ylabel('x')

plt.subplot(1, 3, 2)
plt.title('Displacement in y direction')
plt.imshow(v)
plt.colorbar()
plt.xlabel('y')
plt.ylabel('x')

plt.subplot(1, 3, 3)
plt.title('Displacement in z direction')
plt.imshow(w)
plt.colorbar()
plt.xlabel('y')
plt.ylabel('x')

plt.tight_layout()
plt.savefig('mode_mpl.png', dpi=150)

# Displacement cross sections
plt.figure(2, figsize=(16, 6))
plt.clf()

plt.subplot(1, 3, 1)
plt.plot(x, u[:, int(u.shape[1] / 2)])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('Displacement in x direction')
plt.tight_layout()

plt.subplot(1, 3, 2)
plt.plot(x, v[:, int(v.shape[1] / 2)])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('Displacement in y direction')
plt.tight_layout()

plt.subplot(1, 3, 3)
plt.plot(x, w[:, int(w.shape[1] / 2)])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('Displacement in z direction')
plt.tight_layout()

plt.savefig('cross_section_mpl.png', dpi=150)

# Show plots
plt.show()

I suppose there is something I do not understand when using boxfield with vector data. Can someone help me or show an alternative export of data on a regular grid. I am using fenics 2017.1 (!) in docker under window 10.Thanks in advance.

Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »