Problem when visualising modes with different materials

General
1

Dear all,

I am having a problem when I visualise the modes calculated from linear elasticity theory when I have two different materials. I have followed the tutorial Defining subdomains for different materials — FEniCSx tutorial, where it is defined the DG0 function space for having the given structural properties at the different subdomains.

The only difference with the tutorial is that in my case, my domains are:

Omega0: {x[2] <= 0.6h, x[2] >= 0.6h + h_g}
Omega1: {x[2] >= 0.6h, x[2] <= 0.6h + h_g}

whereas in the tutorial is (with h=1)

Omega0: {x[2] <= 0.5h}
Omega1: {x[2] >= 0.5h}

Here is a minimal code:

import numpy as np
from mpi4py import MPI
import dolfinx
from dolfinx.fem import (form, VectorFunctionSpace, FunctionSpace, Function, dirichletbc, locate_dofs_topological)
from dolfinx import io
from ufl import (inner, div, dx, nabla_grad, sym, Identity, as_tensor, Dn, grad, ds, rhs, lhs, Measure, transpose, as_tensor,
                dot, FacetNormal, indices)
from petsc4py import PETSc
import ufl
from petsc4py.PETSc import ScalarType
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix

gdim = 3
fdim = gdim-1

L, t, h, h_g = 900, 80, 300, 10

Nx = 20
Ny = 10
Nz = 80

mesh = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2*L, t, h])], [Nx, Ny, Nz], mesh.CellType.tetrahedron)

def Omega_0(x):
    return np.logical_or((x[2] <= 0.6*h), (x[2] >= 0.6*h + h_g))

def Omega_1(x):
    return np.logical_and((x[2] >= 0.6*h), (x[2] <= 0.6*h + h_g))

Q = FunctionSpace(mesh, ("DG", 0))
cells_0 = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, Omega_1)

# Material properties
E_0 = 100e9
E_1 = 10e9
nu_0 = 0.23
nu_1 = 0.23
rho_0 = 3000.0
rho_1 = 600.0

# Lame coefficient for constitutive relation
mu_0 = E_0/2./(1+nu_0)
lmbda_0 = E_0*nu_0/(1+nu_0)/(1-2*nu_0)
mu_1 = E_1/2./(1+nu_1)
lmbda_1 = E_1*nu_1/(1+nu_1)/(1-2*nu_1)

mu = Function(Q)
mu.x.array[cells_0] = np.full_like(cells_0, mu_0, dtype=ScalarType)
mu.x.array[cells_1] = np.full_like(cells_1, mu_1, dtype=ScalarType)

rho = Function(Q)
rho.x.array[cells_0] = np.full_like(cells_0, rho_0, dtype=ScalarType)
rho.x.array[cells_1] = np.full_like(cells_1, rho_1, dtype=ScalarType)

lmbda = Function(Q)
lmbda.x.array[cells_0] = np.full_like(cells_0, lmbda_0, dtype=ScalarType)
lmbda.x.array[cells_1] = np.full_like(cells_1, lmbda_1, dtype=ScalarType)

V = VectorFunctionSpace(mesh, ("CG", 1))

def eps(xi):
    return ufl.sym(ufl.grad(xi))

def sigma(xi):
    return 2.0*mu*eps(xi) + lmbda*ufl.nabla_div(xi)*ufl.Identity(len(xi))

xi = ufl.TrialFunction(V)
w = ufl.TestFunction(V)

# Define homogeneous Dirichlet BCs
def left(x):
    return np.isclose(x[0],0)
def right(x):
    return np.isclose(x[0],2*L)
def bottom(x):
    return np.isclose(x[2],0)
def top(x):
    return np.isclose(x[2],h)

left_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, left)
right_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, right)
top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top)
bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom)

u_D = np.array([0, 0, 0], dtype=ScalarType)
bc1 = dirichletbc(u_D, locate_dofs_topological(V, fdim, left_facets), V)
bc2 = dirichletbc(u_D, locate_dofs_topological(V, fdim, right_facets), V)
bc3 = dirichletbc(u_D, locate_dofs_topological(V, fdim, top_facets), V)
bc4 = dirichletbc(u_D, locate_dofs_topological(V, fdim, bottom_facets), V)

bc = [bc1, bc2, bc3, bc4]

# Define stiffness and mass matrices
k_form = form(ufl.inner(sigma(xi),eps(w))*dx)
m_form = form(rho*ufl.inner(xi,w)*dx)

# Assemble matrices
K = assemble_matrix(k_form, bcs=bc, diagonal=628310000)
K.assemble()
M = assemble_matrix(m_form, bcs=bc, diagonal=1/628310000)
M.assemble()

# Create and configure eigenvalue solver
from slepc4py import SLEPc
N_eig = 10
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(N_eig)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
st = SLEPc.ST().create(MPI.COMM_WORLD)
st.setType(SLEPc.ST.Type.SINVERT)
st.setShift(0.1)
st.setFromOptions()
eigensolver.setST(st)
eigensolver.setOperators(K, M)
eigensolver.setFromOptions()

# Compute eigenvalue-eigenvector pairs and save in XDMF file
eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.getVecs()
u_output = Function(V)
u_output.name = "Eigenvector"
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
    with io.XDMFFile(MPI.COMM_WORLD, "eigenvectors.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        for i in range (min(N_eig, evs)):
            l = eigensolver.getEigenpair(i, vr, vi)
            freq = np.sqrt(1000*l.real)/2/np.pi
            print(f"Mode {i}: {freq} Hz")
            u_output.x.array[:] = vr
            xdmf.write_function(u_output, i)

The first mode has this strange shape:

Mode1
Mode11514×795 34.4 KB

In fact, when I define the domains as in the Tutorial, I get a mode that makes sense:

Mode2
Mode21514×795 42.3 KB

Does anyone know what I am doing wrong? It clearly seems that there is something wrong with the definition of the subdomains in my case.

Thanks.

2 
Nz = 80

mesh = mesh.create_box(MPI.COMM_WORLD, [np.array([..., 0]), np.array([..., h])], [..., Nz], ...)

def Omega_0(x):
    return np.logical_or((x[2] <= 0.6*h), (x[2] >= 0.6*h + h_g))

For starters, I would be more careful in the way you define the mesh on the z direction.

Since (Nz = 80)/0.6 is not an integer, you will not have facets exactly at the interface between the two subdomains.

1 Like