Fatal Error in Mesh Refinement

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.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()

b = fe.assemble(L)
delta = fe.PointSource(V, fe.Point(0,0), 1)
fe.solve(A, u.vector(), b)

plt.title("Mesh after refinement")

p = fe.plot(u)

# 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" )

Is it possible that at the 10th refinement iteration a point is added at the origin (or very close to it)? Your expression u_e is singular there.

1 Like

I hadn’t considered that. Thanks