Hello,
I am trying to solve a Stokes equation on a rectangle minus a disk. The mesh was generated like in Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken. I have boundary conditions on the edges of the rectangle and on the boundary of the inner disk and I collect them in the list bcs. When tying to execute the last line, I’m getting Attribute error: ‘list’ object has no attribute ‘_cpp_object’, but I can’t really understand why, because in Stokes equations using Taylor-Hood elements — DOLFINx 0.10.0.0 documentation a very similar thing was done and it worked. I would very much appreciate your help! Sorry for the long code!
from mpi4py import MPI
import pyvista as pv
import numpy as np
from mpi4py import MPI
from dolfinx import default_real_type
from dolfinx.io import XDMFFile
from dolfinx.fem import functionspace, Function, dirichletbc, locate_dofs_geometrical, form, Constant
from dolfinx.fem.petsc import LinearProblem
from ufl import TrialFunction, TestFunction, div, grad, inner, dx
from basix.ufl import element
from petsc4py import PETSc
filename = 'mesh.xdmf'
with XDMFFile(MPI.COMM_WORLD, filename, "r") as xdmf:
msh = xdmf.read_mesh(name="Grid")
P2 = element(
"Lagrange", msh.basix_cell(), degree=2, shape=(msh.geometry.dim,), dtype=default_real_type
)
P1 = element("Lagrange", msh.basix_cell(), degree=1, dtype=default_real_type)
V, Q = functionspace(msh, P2), functionspace(msh, P1)
gdim = 2
def circle_boundary(x):
return np.isclose((x[0] - 10)**2 + (x[1] - 5)**2, 2.5**2)
def left_boundary(x):
return np.isclose(x[0], 0)
def top_bottom_boundary(x):
return np.isclose(x[1], 0) | np.isclose(x[1], 10)
def right_boundary(x):
return np.isclose(x[0], 30)
class InletVelocity():
def __init__(self):
pass
def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = x[1] * (10 - x[1]) / (25.0)
return values
class CircleBoundary():
def __init__(self):
pass
def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = c1 * np.sin(np.arctan2(x[1] - 5, x[0] - 10)) + d1 * np.cos(np.arctan2(x[1] - 5, x[0] - 10)) + c2 * np.sin(np.arctan2(x[1] - 5, x[0] - 10))**2 + d2 * np.cos(np.arctan2(x[1] - 5, x[0] - 10))**2
values[1] = c1 * np.sin(np.arctan2(x[1] - 5, x[0] - 10)) + d1 * np.cos(np.arctan2(x[1] - 5, x[0] - 10)) + c2 * np.sin(np.arctan2(x[1] - 5, x[0] - 10))**2 + d2 * np.cos(np.arctan2(x[1] - 5, x[0] - 10))**2
return values
class NoSlip():
def __init__(self):
pass
def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
return values
class OutletPressure():
def __init__(self):
pass
def __call__(self, x):
values = np.zeros((x.shape[1],), dtype=PETSc.ScalarType)
return values
c1, d1, c2, d2 = 0.1, -0.1, -0.2, 0.2
u_inlet = Function(V)
inlet_velocity = InletVelocity()
u_inlet.interpolate(inlet_velocity)
bcu_inlet = dirichletbc(u_inlet, locate_dofs_geometrical(V, left_boundary))
m = Function(V)
circleboundary = CircleBoundary()
m.interpolate(circleboundary)
bcm = dirichletbc(m, locate_dofs_geometrical(V, circle_boundary))
u_noslip = Function(V)
noslip = NoSlip()
u_noslip.interpolate(noslip)
bcu_noslip = dirichletbc(u_noslip, locate_dofs_geometrical(V, top_bottom_boundary))
p_outlet = Function(Q)
outlet_pressure = OutletPressure()
p_outlet.interpolate(outlet_pressure)
bcp = dirichletbc(p_outlet, locate_dofs_geometrical(Q, right_boundary))
bcs = [bcu_inlet, bcm, bcu_noslip, bcp]
(u, p) = TrialFunction(V), TrialFunction(Q)
(v, q) = TestFunction(V), TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx], [inner(div(u), q) * dx], None])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])
problem = LinearProblem(a, L, [bcu_inlet, bcu_noslip, bcm, bcp])
uh = problem.solve()