Search code examples
pythonnumpymatplotlibscipycomsol

Reflecting 4d symmetric data

I've found out that there is no example of how to reflect symmetric 4d data, which can be very useful when 3d simulations wants to be performed using a symmetric plane to reduce calculations(e.g. ANSYS, COMSOL, etc). This example shows a data file structure corresponding to a COMSOL simulation, which has the structure: X, Y, Z, Amplitude

The model had a symmetry along the Y-plane and was sliced on this plane, so less mesh elements must be calculated. In order to obtain a full Y-plane view (i.e. Y from -0.5 to 0.5), the data has to be reflected along the Y plane.

An example code for such a problem would look like the following:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

### generate a SIMULATION type-file###
X = np.linspace(-0.5, 0.5, 50)
Y = np.linspace(-0.5, 0, 50)    #to be reflected/extended to 0.5
Z = np.linspace(-0.5, 0.5, 50)
Xq, Yq, Zq = np.meshgrid(X, Y, Z)
Amp = 1* np.exp(-((Xq - 0) ** 2 / (0.03) + ( Yq - 0) ** 2 / (0.03) + ( Zq - 0) ** 2 / (0.03)))

datafile = np.vstack([Xq.ravel(), Yq.ravel(), Zq.ravel(), Amp.ravel()]).T   #resemble the simulation data structur, in this case X, Y, Z, Amp

### PYTHON POST-PROCESSING ###
X = datafile[:, 0]
Y = datafile[:, 1]
Z = datafile[:, 2]
Amp = datafile[:, 3]    #Amplitude could be a Pressure, Force, Magnetic field, etc

xq = 0.0    #choose an arbitrary plane to show amplitude distribution over this plane
yq = np.linspace(min(Y), max(Y), 50)  
zq = np.linspace(min(Z), max(Z), 50)
Xq, Yq, Zq = np.meshgrid(xq, yq, zq)
int_plane = griddata((X, Y, Z), Amp, (Xq, Yq, Zq), method='linear')
int_plane = np.squeeze(int_plane)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(Zq, Yq, zs=int_plane)
ax.set_title('3D view');
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

enter image description here An important detail to have in consideration, is that the plane Y=0 must not be duplicated.

The purpose here is to reconstruct the other half of the dataset (i.e. the whole other 3D half space + corresponding amplitudes). A correct reflection of the 3D space would output then a full 3D Gaussian. How could this be efficiently accomplished?

Solution

  • You reflect along the y axis, omitting the last element, and append the reflected data to data itself

    amp = np.vstack((amp, amp[-2::-1]))

    Example:

    In [37]: x, y, z = (-1, 0, 1), (-50, 0), (-5, 0, 5)
    ...: X, Y, X = np.meshgrid(x, y, z)
    ...: data = X+Y+Z
    ...: print("----------------", data, sep='\n')
    ...: data = np.vstack((data, data[-2::-1]))
    ...: print("----------------", data, sep='\n')
----------------
[[[-56 -50 -44]
  [-56 -50 -44]
  [-56 -50 -44]]

 [[ -6   0   6]
  [ -6   0   6]
  [ -6   0   6]]]
----------------
[[[-56 -50 -44]
  [-56 -50 -44]
  [-56 -50 -44]]

 [[ -6   0   6]
  [ -6   0   6]
  [ -6   0   6]]

 [[-56 -50 -44]
  [-56 -50 -44]
  [-56 -50 -44]]]

In [38]: