Yes , it work successfully and just like what I want. Thank you for your help.
And I still need your help.Here are my problem:
- About the boundary defined. When I build boundary without u(y = 0) = 0, the solution was a little different at y = 1 , but totally different at y = 0 , especially the real part. The figure are as follow, first one is without u(y = 0) = 0. The second one is with u(y = 0 ) = 0.The black parts are just because my grid is too dense.
the real part of the numerical solution without u(y = 0) = 0
the imaginary part of the numerical solution without u(y = 0) = 0
the real part of the numerical solution with u(y = 0) = 0
the imaginary part of the numerical solution with u(y = 0) = 0
The difference confuse me. I want you to tell me why. The main code is as follow:
import numpy as np
import math
import ufl
from dolfinx import fem, mesh, plot
from ufl import ds, dx, grad, inner, sin
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (3.0, 1.0)), n=(300, 100),
cell_type=mesh.CellType.triangle,)
V = fem.FunctionSpace(msh, ("Lagrange", 1))
facets = mesh.locate_entities_boundary(msh, dim=1,
marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[0], 3.0))) # without u(y = 0) = 0
# facets = mesh.locate_entities_boundary(msh, dim=1,
# marker=lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0.0),
# np.isclose(x[0], 3.0)),
# np.isclose(x[1], 0.0))) # with u(y = 0) = 0
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
pi = math.pi
x1 = sin(pi * x[0])
y1 = sin(2 * pi * x[1])
f = pi * x1 * (2.5 * y1 + 1j * pi * x[1])
g = (1 + 1j) * sin(pi * x[0])
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx + inner(g, v) * ds
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
And the pyvista part
# uh.real
try:
import pyvista
cells, types, x = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["ur"] = uh.x.array.real
grid.set_active_scalars("ur")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
plotter.show_grid()
plotter.show()
except ModuleNotFoundError:
print("'pyvista' is required to visualise the solution")
print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")
# uh.imag
try:
import pyvista
cells, types, x = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["ui"] = uh.x.array.imag
grid.set_active_scalars("ui")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
plotter.show_grid()
plotter.show()
except ModuleNotFoundError:
print("'pyvista' is required to visualise the solution")
print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")
2.About my recent and initial problem----the answer of “uh * ds”
The answer should be j * 2 / pi , approximately equal to j * 0.636619772
Only with the boundary condition u(y = 0) = 0, I could get
(-2.3239149089002746e-05+0.6365176098632956j)
It’s very close to the analytical solution. That’s sounds a good news for me.
By the way, I want to know how to type formula in the topic and reply.
