Hi everyone,
I’m using FEniCS 2019.1.0 to solve Poisson’s equation on the unit circle with a delta source, i.e.
\nabla^2 u = \delta(\mathbf{x}) \quad \mathrm{on} \; \Omega
\\
u = 0 \quad \mathrm{on} \; \partial\Omega.
Here I’m initially using a coarse mesh and refining it to see if I can get convergence to the analytical solution, u = -\frac{1}{2\pi}\ln(||\mathbf{x}||). However, if I try to refine the mesh more than 10 times (the resulting mesh is not obscenely-fine by my estimation), the program crashes and the following error message is printed to the console,
/home/user/anaconda3/lib/python3.11/site-packages/spyder/plugins/ipythonconsole/scripts/conda-activate.sh: line 18: 18415 Segmentation fault (core dumped) $CONDA_ENV_PYTHON -m spyder_kernels.console -f $SPYDER_KERNEL_SPEC
Fatal Python error: Segmentation fault
Main thread:
Current thread 0x000073d8af23e740 (most recent call first):
File "/home/user/Programming/fenics/mesh_refinement.py", line 50 in <module>
File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/spyder_kernels/py3compat.py", line 356 in compat_exec
File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/spyder_kernels/customize/spydercustomize.py", line 473 in exec_code
File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/spyder_kernels/customize/spydercustomize.py", line 615 in _exec_file
File "/home/user/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/spyder_kernels/customize/spydercustomize.py", line 528 in runfile
File "/tmp/ipykernel_18415/1607777891.py", line 1 in <module>
Restarting kernel...
Find attached the program I’m using. Is there a better way to do mesh refinement that I’m just not aware of?
import fenics as fe
import mshr as msh
import matplotlib.pyplot as plt
N = 5
# (a) using FEniCS built-in 'PointSource' method
pt = fe.Point(0,0)
domain = msh.Circle(pt, 1)
mesh = msh.generate_mesh(domain, N)
plt.figure(1)
fe.plot(mesh)
plt.title("Mesh before refinement")
# Refine mesh at the singularity
cell_markers = fe.MeshFunction("bool", mesh, mesh.topology().dim(), False)
# nor = number of refinements
nor = 25
for n in range(nor):
for cell in fe.cells(mesh):
cell_markers[cell] = cell.contains(pt)
mesh = fe.refine(mesh, cell_markers)
V = fe.FunctionSpace(mesh, "P", 2)
u_D = fe.Expression("0", degree = 2)
def boundary(x, on_boundary):
return on_boundary
bc = fe.DirichletBC(V, u_D, boundary)
u = fe.TrialFunction(V)
v = fe.TestFunction(V)
f = fe.Constant(0)
a = fe.dot( fe.grad(u), fe.grad(v) )*fe.dx
L = f*v*fe.dx
u = fe.Function(V)
assembler = fe.SystemAssembler(a, L, bc)
solver = fe.LUSolver("mumps")
A = fe.Matrix()
solver.set_operator(A)
assembler.assemble(A)
b = fe.assemble(L)
delta = fe.PointSource(V, fe.Point(0,0), 1)
delta.apply(b)
bc.apply(b)
fe.solve(A, u.vector(), b)
plt.figure(2)
fe.plot(mesh)
plt.title("Mesh after refinement")
plt.figure()
plt.title("PointSource")
p = fe.plot(u)
p.set_cmap("viridis")
plt.colorbar(p)
# Compute error
u_e = fe.Expression('-1/(2*3.14159)*'
+ 'std::log( pow( pow(x[0],2) + pow(x[1],2), 1/2) )',
degree = 6)
V = u.function_space()
mesh = V.mesh()
degree = V.ufl_element().degree()
W = fe.FunctionSpace(mesh, "P", degree + 3)
u_e_w = fe.interpolate(u_e, W)
u_w = fe.interpolate(u, W)
e_w = fe.Function(W)
e_w.vector()[:] = u_e_w.vector().get_local() - u_w.vector().get_local()
error = e_w**2*fe.dx
E = fe.errornorm(u_e, u, 'L2')
print("error in exact delta = ", E, "\n" )