Refining Mesh Close to a Point

mesh
1

Hi everyone,

I’m using FEniCS 2019.1.0 to refine a mesh successively in the neighborhood of some point in the domain. What I do in the code snippet below is mark each cell for refinement having at least one vertex a distance less than eps away from the point of interest.

When I try to do this, however, my program crashes, even for rather modest values of Num_refinements. Is there a way to salvage this?

from fenics import *
import numpy as np

pt = np.array([0.5, 0.5])


# Create mesh 
nx = ny = 2
mesh = UnitSquareMesh(nx, ny)

Num_refinements = 10

refine_cell = MeshFunction("bool", mesh, mesh.topology().dim())

eps = 0.01
for i in range(Num_refinements):
    ctr = 0
    for c in cells(mesh):
        ctr = ctr+1
        dists = np.zeros(3)
        conds = np.zeros(3)
        vertices_x = sorted(c.get_vertex_coordinates()[::2])
        vertices_y = sorted(c.get_vertex_coordinates()[1::2])
        for i in range(3):
            dists[i] = np.sqrt( (vertices_x[i] - pt[0])**2\
                               + (vertices_y[i] - pt[1])**2 )
            conds[i] = dists[i] < eps
        if any(conds):
            refine_cell[c] = True
        else:
            refine_cell[c] = False
    mesh = refine(mesh, refine_cell)
2

refine_cell should be inside the refinement loop, as you keep overwriting the mesh variable:

from fenics import *
import numpy as np

pt = np.array([0.5, 0.5])


# Create mesh 
nx = ny = 2
mesh = UnitSquareMesh(nx, ny)

Num_refinements = 10


eps = 0.01
for i in range(Num_refinements):
    refine_cell = MeshFunction("bool", mesh, mesh.topology().dim())
    ctr = 0
    for c in cells(mesh):
        ctr = ctr+1
        dists = np.zeros(3)
        conds = np.zeros(3)
        vertices_x = sorted(c.get_vertex_coordinates()[::2])
        vertices_y = sorted(c.get_vertex_coordinates()[1::2])
        for i in range(3):
            dists[i] = np.sqrt( (vertices_x[i] - pt[0])**2\
                               + (vertices_y[i] - pt[1])**2 )
            conds[i] = dists[i] < eps
        if any(conds):
            refine_cell[c] = True
        else:
            refine_cell[c] = False
    mesh = refine(mesh, refine_cell)
1 Like