Hi,
I try to solve a simple 3D Navier-Stokes 3D (sphere) problem.
And for this, I have considered this example:Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial
However, I am more used to FEniCS 2019, where I just mark the BC. But, in FEniCS-X ns_code2 example, the BC is marked with facets.
I have managed to read the mesh and create the function space in FEniCS-X
import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx.cpp.mesh import CellType
from dolfinx.io import (XDMFFile, extract_gmsh_geometry,
extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh)
from dolfinx.fem import (Constant, DirichletBC, Function, FunctionSpace, apply_lifting, assemble_matrix,
assemble_scalar, assemble_vector, create_vector, locate_dofs_topological, set_bc)
from petsc4py import PETSc
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
print(f"DOLFINx version: {dolfinx.__version__} is installed")
filename="sphere-2.xdmf"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "sphere-2.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
t = 0
T = 1 #8 # Final time
dt = 1/1600 # Time step size
num_steps = int(T/dt)
k = Constant(mesh, dt)
mu = Constant(mesh, 0.001) # Dynamic viscosity
rho = Constant(mesh, 1) # Density
v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)
fdim = mesh.topology.dim
For a similar problem in FEniCS 2019, I have set the BC:
# Define function spaces for pressure and velocity
V = VectorFunctionSpace(mesh, "Lagrange", degree=2, dim=3)
Q = FunctionSpace(mesh, "Lagrange", 1)
class Inflow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[0] - 0.) < DOLFIN_EPS)
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[0] - 10.) < DOLFIN_EPS)
class Walls(SubDomain):
def inside(self, x, on_boundary):
result = on_boundary and ((abs(x[1] + 0.) < DOLFIN_EPS) or (abs(x[1] - 5.) < DOLFIN_EPS) or
(abs(x[2] - 0.) < DOLFIN_EPS) or (abs(x[2] - 5.) < DOLFIN_EPS))
print(result)
return result
class Sphere(SubDomain):
def inside(self, x, on_boundary):
result = on_boundary\
and x[0] > 2.0 - DOLFIN_EPS and x[0] < 4.0 + DOLFIN_EPS\
and x[1] > 2.0 - DOLFIN_EPS and x[1] < 4.0 + DOLFIN_EPS\
and x[2] > 2.0 - DOLFIN_EPS and x[2] < 4.0 + DOLFIN_EPS
return result
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() , 0)
Inflow().mark(boundaries, 1)
Outflow().mark(boundaries, 2)
Walls().mark(boundaries, 3)
Sphere().mark(boundaries, 4)
inflow_profile = ('1.', '0.', '0.')
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=3), boundaries, 1)
bcp_outflow = DirichletBC(Q, Constant(0.), boundaries, 2)
bcu_walls = DirichletBC(V, Constant((0., 0., 0.)), boundaries, 3)
bcu_sphere = DirichletBC(V, Constant((0., 0., 0.)), boundaries, 4)
bcu = [bcu_inflow, bcu_walls, bcu_sphere]
bcp = [bcp_outflow]
Could you please tell me how to do it in a similar way in FEniCS-X or is there any other smart way to do set the BC?
