You should not apply your bcs to the collapsed subspaces V
and X
, see:
In your case, changing your first MWE to the following code should allow you to solve using MixedElement
s. Note that I corrected the definitions of right_boundary
and p0
to match the definitions in your second example.
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import FunctionSpace, Function, dirichletbc, locate_dofs_topological, Constant
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType
import ufl
import numpy as np
def bottom_boundary(x):
return np.isclose(x[1],0)
def left_boundary(x):
return np.isclose(x[0],0)
def right_boundary(x):
return np.isclose(x[0],1)
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u,Lame_params):
return Lame_params[0]*ufl.nabla_div(u)*ufl.Identity(u.geometric_dimension()) + 2*Lame_params[1]*epsilon(u)
#Main simuluation
rho = 2000 #density
Ym = 2e7 #Young modulus
Pr = 0.3 #Poisson ratio
g = 9.8 #accel due to gravity
rho_w = 1000 #density of water
Lame_params = (Ym*Pr)/((1+Pr)*(1-2*Pr)), 0.5*Ym/(1+Pr)
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.quadrilateral)
V_el = ufl.VectorElement("CG", domain.ufl_cell(), 1)
X_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
M_el = ufl.MixedElement([V_el, X_el])
VX = FunctionSpace(domain,M_el)
(u,p) = ufl.TrialFunctions(VX)
(v,q) = ufl.TestFunctions(VX)
bcs = []
#boundary conditions on displacement
V, V_dofs = VX.sub(0).collapse()
X, X_dofs = VX.sub(1).collapse()
fdim = domain.topology.dim - 1
p0 = Function(X)
p0.interpolate(lambda x: -rho_w*g*x[1])
u0 = Function(V)
boundary_facets = mesh.locate_entities_boundary(domain, fdim, bottom_boundary)
bcs.append(dirichletbc(u0,locate_dofs_topological((VX.sub(0), V), fdim, boundary_facets), VX.sub(0)))
bcs.append(dirichletbc(p0,locate_dofs_topological((VX.sub(1), X), fdim, boundary_facets), VX.sub(1)))
boundary_facets = mesh.locate_entities_boundary(domain, fdim, left_boundary)
bcs.append(dirichletbc(u0,locate_dofs_topological((VX.sub(0), V), fdim, boundary_facets), VX.sub(0)))
bcs.append(dirichletbc(p0,locate_dofs_topological((VX.sub(1), X), fdim, boundary_facets), VX.sub(1)))
boundary_facets = mesh.locate_entities_boundary(domain, fdim, right_boundary)
bcs.append(dirichletbc(p0,locate_dofs_topological((VX.sub(1), X), fdim, boundary_facets), VX.sub(1)))
G = Constant(domain,ScalarType((0,-g)))
a = ufl.inner(sigma(u,Lame_params),epsilon(v))*ufl.dx
a += -ufl.inner(p*ufl.Identity(u.geometric_dimension()),epsilon(v))*ufl.dx
a += ufl.dot(ufl.grad(p),ufl.grad(q))*ufl.dx
L = ufl.dot(rho*G,v)*ufl.dx + Constant(domain,ScalarType(0))*q*ufl.dx
problem = LinearProblem(a,L,bcs,petsc_options={"ksp_type":"preonly","pc_type":"lu"})
up = problem.solve()
u,p = up.split()
print(u.x.array)
print(p.x.array)
from dolfinx import io
with io.XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(u)
with io.XDMFFile(MPI.COMM_WORLD, "p.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(p)