Choose Refinement of Cell and Interpolation of Functions on New Mesh

mesh
1

Fenics Version: 2019.2.0.dev0

Hello, I am trying to refine a 1D spatial mesh cell into n new 1D spatial cells. I want to take advantage of known asymptotic error decay to reduce the number of times the mesh needs to be refined to hit some local error tolerance. When I use the refine() function, it refines the 1D spatial cell using bisection and splits the cell in two. I want to be able to split the cell into however many equal cells that I want. I was wondering if there was any inbuilt way in Fenics to do this.

On a separate note, I was wondering if there are any issues interpolating function from a coarse spatial mesh to a commonly refined spatial mesh in Fenics. When I look at the mesh pre- and post-refinement I see some differences that concern me, but the output of specific numerical values appears the same. Specifically, the order of the cells seems fine, but the order of the coordinates of the vertices seems strange with the midpoint of the two new cells at the end rather than in order. Also, the functions when evaluated at common points appear to give the same output. But the graphs of the function after mesh refinement is strange, it appears to be graphing backwards towards the refined cell. So I’m not sure if there are any potential issues, but that’s why I’m asking.

My minimum working example of what I am seeing.

from fenics import *
import dolfin as df
set_log_level(LogLevel.WARNING)

Nx = 5
deg = 1
mesh = UnitIntervalMesh(Nx)
V = FunctionSpace(mesh, 'CG', deg)
f = Expression('sin(pi*x[0])', degree = 6)
f_func = interpolate(f, V)
plot(mesh)
mesh.cells()
mesh.coordinates()
plot(f_func)

cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
i = 0
for cell in cells(mesh):
	if i == 0: #Placeholder criterion
		cell_markers[cell] = True
	i = i + 1

mesh_new = refine(mesh, cell_markers)
plot(mesh_new)
mesh_new.cells()
mesh_new.coordinates()

V_new = FunctionSpace(mesh_new, 'CG', deg1)
f_func_new = interpolate(f_func, V_new)
plot(f_func_new)

f_func(0.0)
f_func(0.1)
f_func(0.2)
f_func(0.3)
f_func(0.4)
f_func(0.5)
f_func(0.6)
f_func(0.7)
f_func(0.8)
f_func(0.9)
f_func(1.0)

f_func_new(0.0)
f_func_new(0.1)
f_func_new(0.2)
f_func_new(0.3)
f_func_new(0.4)
f_func_new(0.5)
f_func_new(0.6)
f_func_new(0.7)
f_func_new(0.8)
f_func_new(0.9)
f_func_new(1.0)
2

There is no built-in way of doing that (except repeated refinement in a given area)

A mesh consists of two things, the

  1. geometry (mesh.coordinates()) which is a list of the physical coordinates.
  2. topology. The connectivity from any cell to a set of vertices (mesh.cells()), as you can see, by calling
for cell in cells(mesh_new):
    print(mesh_new.coordinates()[mesh_new.cells()[cell.index()]])

you get the vertices of each cell:

[[0. ]
 [0.1]]
[[0.2]
 [0.1]]
[[0.2]
 [0.4]]
[[0.4]
 [0.6]]
[[0.6]
 [0.8]]
[[0.8]
 [1. ]]

which shows that they are correctly linked.

In order, as in from lowest to highest value, is not an efficient way of refining the mesh (due to memory allocation). Thus it is alot faster to add new vertices at the end than in the middle of an array.

This is just due to how the built-in plot module plots 1D meshes.
A better way for visualizing is by sorting the dofs when plotting

vals = f_func_new.compute_vertex_values()
coords = mesh_new.coordinates()
sort_order = np.argsort(coords[:, 0])
plt.plot(coords[sort_order, 0], vals[sort_order], "-ro")

resulting in
new_func

1 Like
3

Thank you for the response. It was very helpful and now I have some new questions.

  1. To make sure I understand, when I refine a mesh, even though the list of mesh nodes is not “in order”, overall it is still ordered as the connections or topology of the nodes are tracked?

  2. Does this means that the cells are ordered? From my testing this appears so since I can use the following code added to my previous code and change the criterion for refinement to select any cell in order to be refined. In other words, i == 0 selects the first cell in the graphed mesh to be refined, i == 1 is the second cell, etc. up to i == 5 for the last cell. Does this process hold in general?

i = 0
for cell in cells(mesh_new):
	if i == 6: #Placeholder criterion
		cell_markers_new[cell] = True
	i = i + 1

mesh_new_2 = refine(mesh_new, cell_markers_new)
plot(mesh_new_2)
  1. I saw that there was an adapt() function on other threads I searched. I don’t know if it is necessary to use the adapt() function vs the interpolate() function when transitioning f_func to f_func_new. Or what the differences are. From the documentation, it looks like it behaves like refine() when passed a mesh. Is it okay to just do the following on meshes that are common refinements of the previous mesh? Or are there issues that might crop up in Fenics?
f_func_new = interpolate(f_func, mesh_new)
4

Yes

This does not hold in general. There is no special implementation of a 1D unit interval in Dolfin when it comes to refinement. The mesh and refinement procedure is the same on any input mesh, thus it has to hold in general, as the ordering you talk about only makes sense on intervals

1 Like
5

Awesome! Thank you for helping me with my mesh questions.

My only question left is:

Should I use interpolate() to move f_func from the old mesh, mesh, to the new refined mesh or should I use a function like adapt()? I’ve seen a bunch of old threads that use adapt(), but comparing the old, 1.5.0 docs, to the 2019.1.0 docs and it looks like a bunch of adapt() functionality changed/was removed. Using interpolate() in my code appears to work fine, but I would like to know if this is only true for my specific example. Is there any case where I would want to use a function other than interpolate() to move f_func from mesh to mesh_new?

6

Adapt is used for meshfunctions, not functions.

You should use interpolate to move functions between meshes.

1 Like
7

Thank you for all the help! Marking as solved!