Relevant DOLFINx modules:
dolfinx.mesh
: Classes and functions related to the computational domaindolfinx.fem
: Finite element method functionalitydolfinx.io
: Input/Output (read/write) functionalitydolfinx.plot
: Convenience functions for exporting plotting datadolfinx.la
: Functions related to linear algebra structures (matrices/vectors)from dolfinx import mesh, fem, io, plot, la
from mpi4py import MPI
length, height = 10, 3
Nx, Ny = 80, 60
extent = [[0., 0.], [length, height]]
domain = mesh.create_rectangle(
MPI.COMM_WORLD, extent, [Nx, Ny], mesh.CellType.quadrilateral)
local_domain = mesh.create_rectangle(
MPI.COMM_SELF, extent, [Nx, Ny], mesh.CellType.quadrilateral)
with $u_D = y\cos(0.25t)$, $f=0$. For this example, we take $\Omega$ to be rectangle defined above, $\Omega_\text{D}$ if the left-hand edge of the rectangle, and $\Omega_\text{N}$ is the remaining three edges of the rectangle.
from ufl import (TestFunction, SpatialCoordinate, TrialFunction,
as_vector, dx, grad, inner, system)
V = fem.FunctionSpace(domain, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
un = fem.Function(V)
f = fem.Constant(domain, 0.0)
mu = fem.Constant(domain, 2.3)
dt = fem.Constant(domain, 0.05)
F = inner(u - un, v) * dx + dt * mu * inner(grad(u), grad(v)) * dx
F -= dt * inner(f, v) * dx
(a, L) = system(F)
import numpy as np
def uD_function(t):
return lambda x: x[1] * np.cos(0.25 * t)
uD = fem.Function(V)
t = 0
uD.interpolate(uD_function(t))
def dirichlet_facets(x):
return np.isclose(x[0], length)
tdim = domain.topology.dim
bc_facets = mesh.locate_entities_boundary(
domain, tdim - 1, dirichlet_facets)
bndry_dofs = fem.locate_dofs_topological(V, tdim - 1, bc_facets)
bcs = [fem.dirichletbc(uD, bndry_dofs)]
compiled_a = fem.form(a)
A = fem.petsc.assemble_matrix(compiled_a, bcs=bcs)
A.assemble()
compiled_L = fem.form(L)
b = fem.Function(V)
from petsc4py import PETSc
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.CG)
pc = solver.getPC()
pc.setType(PETSc.PC.Type.HYPRE)
pc.setHYPREType("boomeramg")
import pyvista
import matplotlib.pyplot as plt
pyvista.start_xvfb(0.5) # Start virtual framebuffer for plotting
plotter = pyvista.Plotter()
plotter.open_gif("u_time.gif")
topology, cells, geometry = plot.create_vtk_mesh(V)
uh = fem.Function(V)
grid = pyvista.UnstructuredGrid(topology, cells, geometry)
grid.point_data["uh"] = uh.x.array
renderer = plotter.add_mesh(grid, show_edges=True, lighting=False,
cmap=viridis, scalar_bar_args=sargs,
clim=[0, height])
T = np.pi
while t < T:
# Update boundary condition
t += dt.value
uD.interpolate(uD_function(t))
# Assemble RHS
b.x.array[:] = 0
fem.petsc.assemble_vector(b.vector, compiled_L)
# Apply boundary condition
fem.petsc.apply_lifting(b.vector, [compiled_a], [bcs])
b.x.scatter_reverse(la.ScatterMode.add)
fem.petsc.set_bc(b.vector, bcs)
# Solve linear problem
solver.solve(b.vector, uh.vector)
uh.x.scatter_forward()
# Update un
un.x.array[:] = uh.x.array
# Update plotter
plotter.update_scalars(uh.x.array, render=False)
plotter.write_frame()
x = SpatialCoordinate(domain)
x_grad = inner(as_vector((x[1], x[0])), grad(uh))
W = fem.FunctionSpace(domain, ("DQ", 1))
expr = fem.Expression(x_grad, W.element.interpolation_points())
w = fem.Function(W)
w.interpolate(expr)