Getting [inf inf ... inf] in this coupled displacement-pressure problem implemented with MixedElements

variational formulation
2

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 MixedElements. 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)
3 Likes