Refining mesh with marked edges using locate_entities

[alex24]

Hi,

I am trying to refine a 3D mesh by checking which edges are within a sphere with some user defined center coordinates and radius. This code is based on an earlier code using MeshFunction:

Local Mesh Refinement

Hi @zaccandels, you can also pass a MeshFunction to refine which indicates whether a given cell should be refined or not. This then needs to local refinement and you just need to write some code for a heuristic that refines close to your interior point. Here is an example of what this could like, but I decided to refine all triangles that interest the vertical lines x = 2.5 and x = 3.5. refine_cell = MeshFunction("bool", mesh, mesh.topology().dim()) for c in cells(mesh): vertices_x = sorted…

For some reason the code crashes when I try refining the mesh unifomly twice before the local refinement (num_refinements is set to 2), and using VaraView to view the local refinement does not appear to show the refinement happening in the right cells. Am I using the function location_to_mark() correctly within mesh.locate_entities()? Is the function finding all the edges in the volume of the mesh? Are the x variables the 2 nodes of each edge?

Thank you,
Alex

The code is shown below:

from mpi4py import MPI
from dolfinx import mesh, io
import basix

import numpy as np


# Simple cube mesh
points = np.array([[0.0, 0.0, 0.0],
                   [0.0, 1.0, 0.0],
                   [0.0, 1.0, 1.0],
                   [0.0, 0.0, 1.0],
                   [1.0, 0.0, 0.0],
                   [1.0, 1.0, 0.0],
                   [1.0, 1.0, 1.0],
                   [1.0, 0.0, 1.0]])
cells  = np.array([[0, 1, 3, 4],
                   [7, 3, 4 ,6],
                   [1, 2, 3, 6],
                   [3, 1, 6, 4],
                   [5, 1, 6, 4]])

finiteElement = basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3, ))

domain = mesh.create_mesh(MPI.COMM_WORLD, cells, points, finiteElement)




fdim = domain.topology.dim - 1


num_refinements = 2
for n in np.arange(num_refinements):
    domain.topology.create_entities(1)
    domain = mesh.refine(domain)
    domain.topology.create_connectivity(fdim, fdim+1)



s_c = np.array([0.0, 0.0, 0.0])

refine_radius = 1.0


def location_to_mark(x):
    
    return np.sqrt((x[0] - s_c[0])**2 + (x[1] - s_c[1])**2 + (x[2] - s_c[2])**2) < (refine_radius + 1.0e-6)


marked_edges = mesh.locate_entities(domain, fdim, location_to_mark)

domain.topology.create_entities(1)
domain = mesh.refine(domain, marked_edges)

domain.topology.create_connectivity(fdim, fdim+1)


xdmf_vec= io.XDMFFile(domain.comm, 'mesh.xdmf', "w")
xdmf_vec.write_mesh(domain)
xdmf_vec.close()

[dokken]

@alex24

marked_edges = mesh.locate_entities(domain, fdim, location_to_mark)

In 3D these are not edges, but facets.
Change this to

marked_edges = mesh.locate_entities(domain, 1, location_to_mark)

[alex24]

I got it thank you; the 1 is the dimension of an edge