Reordering nodal values in accordance with meshgrid created using numpy

dolfinx
,
1

After extracting nodal values from the solution array, I want to reorder these values so that they correspond to new_coordinates. If I use the numpy array instead of dof_coor it is working well, I am not sure what is going wrong if I use dof_coor array to search for coordinates. It is giving key missing error. I would appreciate any help.

import numpy as np
import matplotlib.pyplot as plt
from dolfinx.fem import functionspace
from dolfinx.mesh import create_unit_square
from mpi4py import MPI

# define mesh and function space
N = 3
degree = 1
msh = create_unit_square(MPI.COMM_WORLD,N,N)
V = functionspace(msh, ("Lagrange", degree))
x = SpatialCoordinate(msh)
uf = fem.Function(V)
uf.interpolate(lambda x : x[0] + x[1])

# points = msh.geometry.x # I can use either this or tabulate_dof_coordinates() to get coordinates of corresponding nodal values
# dof_coor = np.delete(points,2,1)
array = uf.x.array
print(array)
# print(dof_coor.dtype)
dof_coor = V.tabulate_dof_coordinates()
dof_coor = np.delete(dof_coor,2,1)
print(dof_coor)

import numpy as np
import matplotlib.pyplot as plt

# Original coordinates and values if I use this with new_coords then it works well
# old_coords = np.array([[0.66666667, 0.        ],
#  [1.        , 0.        ],
#  [1.        , 0.33333333],
#  [0.66666667, 0.33333333],
#  [0.33333333, 0.        ],
#  [1.        , 0.66666667],
#  [0.33333333, 0.33333333],
#  [0.66666667, 0.66666667],
#  [0.        , 0.        ],
#  [1.        , 1.        ],
#  [0.        , 0.33333333],
#  [0.33333333, 0.66666667],
#  [0.66666667, 1.        ],
#  [0.        , 0.66666667],
#  [0.33333333, 1.        ],
#  [0.        , 1.        ]])
old_coords = dof_coor # but if I use this and meshgrid coordinates it doesn't work well

values = array

# Given new coordinates  
# new_coords = np.array([[0.        , 0.        ],
#                        [0.33333333, 0.        ],
#                        [0.66666667, 0.        ],
#                        [1.        , 0.        ],
#                        [0.       , 0.33333333],
#                        [0.33333333, 0.33333333],
#                        [0.66666667, 0.33333333],
#                        [1.        , 0.33333333],
#                        [0.        , 0.66666667],
#                        [0.33333333, 0.66666667],
#                        [0.66666667, 0.66666667],
#                        [1.        , 0.66666667],
#                        [0.        , 1.        ],
#                        [0.33333333, 1.        ],
#                        [0.66666667, 1.        ],
#                        [1.        , 1.        ]])
# mesh
a = 0
b = 1

x, y = np.linspace(a,b,4,dtype = np.float64), np.linspace(a,b,4,dtype = np.float64)
X,Y = np.meshgrid(x,y)
X_flat, Y_flat = X.flatten(), Y.flatten()
new_coords = np.vstack((X_flat, Y_flat)).T
print(new_coords)
# Create a mapping from old coordinates to values
mapping = {tuple(coord): value for coord, value in zip(old_coords, values)}
print(mapping)
# Use the mapping to assign values to new coordinates
new_values = np.array([mapping[tuple(coord)] for coord in new_coords])
new_values = new_values.reshape((4,4))
X,Y = np.meshgrid(new_values[:,0],new_values[:,1])
plt.scatter(X,Y)

error

KeyError                                  Traceback (most recent call last)
<ipython-input-92-ef8cf7fe0e24> in <cell line: 56>()
     54 print(mapping)
     55 # Use the mapping to assign values to new coordinates
---> 56 new_values = np.array([mapping[tuple(coord)] for coord in new_coords])
     57 new_values = new_values.reshape((4,4))
     58 X,Y = np.meshgrid(new_values[:,0],new_values[:,1])

<ipython-input-92-ef8cf7fe0e24> in <listcomp>(.0)
     54 print(mapping)
     55 # Use the mapping to assign values to new coordinates
---> 56 new_values = np.array([mapping[tuple(coord)] for coord in new_coords])
     57 new_values = new_values.reshape((4,4))
     58 X,Y = np.meshgrid(new_values[:,0],new_values[:,1])

KeyError: (0.0, 0.0)
2

Hi Aijaz,

Your problem here lies with the representation of zero for floating numbers. I ran your code and inspection of old_coords shows that (on my computer) the point (0.0, 0.0) was represented as (7.18076024e-18, -5.55608928e-18), so the operation mapping[tuple(old_coords)] won’t work because mapping has the key (7.18076024e-18, -5.55608928e-18) and not (0.0, 0.0). So you would have to come up with another way of remapping the coordinates. From the looks of your new_coords it seems you could e.g. try to sort your list of values based on the y-coordinates of their coordinates.

Just as a sidenote and a comment to your comment in the code that you can use either msh.geometry.x or tabulate_dof_coordinates() to get the nodal values: the fact that the shape of tabulate_dof_coordinates() (the finite element dof coordinates/the nodal values) matches the shape of msh.geometry.x (the mesh coordinates) in your example is due to your elements being Lagrange polynomials of order 1, in which case the dofs coincide with your mesh triangularization. If you would e.g. use order 2, you would have more dof coordinates than mesh coordinates, and in turn tabulate_dof_coordinates() and msh.geometry.x would not be the same. The finite element nodal values would then be in tabulate_dof_coordinates(), but not in msh.geometry.x.

Hope this is of help :slight_smile:

Best,
Halvor

2 Likes
3

Consider the following minimal example:

import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
import dolfinx
import ufl

# define mesh and function space
N = 10
degree = 1
msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD,N,N)
V = dolfinx.fem.functionspace(msh, ("Lagrange", degree))
x = ufl.SpatialCoordinate(msh)
uf = dolfinx.fem.Function(V)
uf.interpolate(lambda x : np.sin(2*np.pi*x[0]) + 3*x[1])

# points = msh.geometry.x # I can use either this or tabulate_dof_coordinates() to get coordinates of corresponding nodal values
# dof_coor = np.delete(points,2,1)
array = uf.x.array
# print(dof_coor.dtype)
dof_coord = V.tabulate_dof_coordinates()[:, :2]

# Compute input index
indices = dof_coord / (1./N)
input_indices = np.round(indices).astype(np.int32)
y_sort = np.argsort(input_indices[:, 1])
sorted_coords = dof_coord[y_sort]
sorted_values = array[y_sort]
# Check that x-axis is sorted by increasing x coordinate (per y coord)
for i in range(N+1):
    sub_indices = np.flatnonzero(input_indices[y_sort, 1]==i)
    sub_sort = np.argsort(input_indices[y_sort[sub_indices], 0])
    sub_coords = sorted_coords[sub_indices]
    sub_vals = sorted_values[sub_indices]
    sorted_coords[sub_indices] = sub_coords[sub_sort]
    sorted_values[sub_indices] = sub_vals[sub_sort]



sorted_values = sorted_values.reshape(N+1, N+1)
X_ = sorted_coords[:, 0].reshape(N+1, N+1)
Y_ = sorted_coords[:, 1].reshape(N+1, N+1)

a = 0
b = 1

x, y = np.linspace(a,b,N+1,dtype = np.float64), np.linspace(a,b,N+1,dtype = np.float64)
X,Y = np.meshgrid(x,y)


np.testing.assert_allclose(X_, X, atol=1e-15, rtol=1e-12)
np.testing.assert_allclose(Y_, Y, atol=1e-15, rtol=1e-12)

X_flat, Y_flat = X.flatten(), Y.flatten()
# Use the mapping to assign values to new coordinates
fig = plt.figure()
ax = plt.gca()

im1 = ax.imshow(sorted_values, cmap="viridis")
ax.invert_yaxis()
fig.colorbar(im1)
plt.savefig("test.png")

fig2 = plt.figure()
ax2 = plt.gca()
im2 = ax2.scatter(dof_coord[:,0], dof_coord[:,1], c=uf.x.array, cmap="viridis")

fig2.colorbar(im2)
plt.savefig("test2.png")

yielding:
test
test2

2 Likes
4

@hherlyng Thanks for your comments.

5

@dokken Thank you, appreciated as always!