Logarithmically spaced mesh nodes

mesh
1

Probably a straightforward question, but a quick Google search didn’t yield any results so I’m asking it here.

How would one create a rectangular mesh with one of the axis logarithmically spaced?
For example, I want the nodes on the x-axis to be separated linearly with, e.g., np.linspace(0, 100, 21); meanwhile, the y-axis nodes to be at locations that are separated as if they are on a log-axis np.logspace(np.log10(1), np.log10(1000), 21)

2

Hi, I made the custom mesh using Expression() like this -

from dolfin import *
import numpy as np

# Define the number of intervals and the range of the x-axis
num_intervals_x = 20
x_start = 0.0
x_end = 100.0

# Define the number of intervals and the range of the y-axis
num_intervals_y = 20
y_start = 1.0
y_end = 1000.0

# Create the linearly spaced mesh for the x-axis
mesh_x = IntervalMesh(num_intervals_x, x_start, x_end)

# Create the logarithmically spaced mesh for the y-axis
log_space = Expression("exp(log_start + (log_end - log_start)*x[0]/L)", degree=1,
                       log_start=np.log(y_start), log_end=np.log(y_end), L=x_end)

# Function space on the x-axis mesh
V = FunctionSpace(mesh_x, "CG", 1) 
mesh_y = interpolate(log_space, V)

# Get the coordinates of the mesh for the x-axis
coords_x = mesh_x.coordinates()

# Create the rectangular mesh using the product of the x and y meshes
mesh = RectangleMesh(Point(coords_x.min(), y_start),
                     Point(coords_x.max(), y_end),
                     num_intervals_x, num_intervals_y)

You can also verify the coordinates of the mesh -

print(mesh.coordinates())
3

@violetus Thanks, but the code still seems to create a mesh with a linearly spaced y-axis. The output of the y mesh coordinates is:

>>> print(mesh.coordinates()[:, 1][::21])
[   1.     50.95  100.9   150.85  200.8   250.75  300.7   350.65  400.6
  450.55  500.5   550.45  600.4   650.35  700.3   750.25  800.2   850.15
  900.1   950.05 1000.  ]

Would it be wrong to just edit the mesh coordinates after creating them with something as simple as:

mesh = RectangleMesh(Point(0, 0), Point(100, 4), 20, 20)
for i_point in range(len(mesh.coordinates())):
    mesh.coordinates()[i_point] = [mesh.coordinates()[i_point][0], 
                                   10**mesh.coordinates()[i_point][1]]
4

Hi @DrKorvin, thanks for correcting me, I misunderstood the question and thought that what I did is correct. I ran your code and visualized it -

image
image1917×865 161 KB

Looks fine to me.

5

@violetus Sounds good. Many thanks for the help!

1 Like