IntuitionRF - Mesh your OpenEMS Simulations with Blender !

[Juleinn]

Hello !

I recently saw the Rotman Lens tutorial by @samerps and I though I wanted to do the same :)
Now meshing the simulation correctly can be quite time consuming, especially when its comes to debugging PEC dumps of STL files or arbitrary shapes.
So I made a small plugin for Blender that lets you define the meshlines from within blender, preview the smoothmeshlines, and debug the PEC (overlap it on the existing antenna/trace design).

I hope you can find usefull !

Here is the plugin : https://github.com/Juleinn/IntuitionRF

And there is a demo/tutorial on Youtube : https://www.youtube.com/watch?v=oCE_hrCGen4

It shows a example of a Patch Antenna, a PIFA antenna, and a trace on a PCB, made with KiCAD

[samerps]

excellent work, thank you very much for creating this, I will definitely experiment with it!

[Juleinn]

I'll make sure to look into that !

The surface currents slices don't show anything, even after scaling. I put the slice both exactly on the PEC layer and above/below with no success.

Sure, you need to inherit bpy.types.Node to create the node. You can then add sockets to the node, either existing ones (geometry) or create custom socket types by extending bpy.types.NodeSocket. Registering the node alone won't let you add it to a node graph, you'll need to add it to the menu like this :

import bpy
from bpy.types import Node, NodeSocket, Operator
from bpy.props import PointerProperty


# Custom socket type
class RFDomainSocket(NodeSocket):
    bl_idname = 'RFDomainSocketType'
    bl_label = 'RF Domain'

    def draw(self, context, layout, node, text):
        layout.label(text=text)

    def draw_color(self, context, node):
        return (1.0, 0.4, 0.216, 0.5)  # RGBA
    
# Custom socket type for Object input
class MyObjectSocket(NodeSocket):
    bl_idname = 'MyObjectSocketType'
    bl_label = 'Object Socket'

    def draw(self, context, layout, node, text):
        layout.prop(self, "object", text=text)

    def draw_color(self, context, node):
        return (0.6, 0.6, 0.2, 1.0)  # RGBA

    @classmethod
    def from_object(cls, node, object):
        return cls(node=node, object=object)

    # Ensure the socket type accepts Object type
    object: PointerProperty(name="Object", type=bpy.types.Object)

# Custom node type
class RFSimNode(Node):
    bl_idname = 'RFSimNodeType'
    bl_label = 'RF Sim'

    def init(self, context):
        self.inputs.new('RFDomainSocketType', "Meshing lines")

    def update(self):
        pass

    def draw_buttons(self, context, layout):
        layout.label(text="Run SIM")
        layout.operator("node.runsim_button", text="Run SIM")
        layout.operator("node.previewpec_button", text="Preview PEC")

    def draw_label(self):
        return "RF Simulation"
    
class RFSimVertexGroupToMeshLinesNode(Node):
    bl_idname = 'RFSimVertexGroupToMeshLinesNodeType'
    bl_label = 'Vertex Group to mesh lines'

    def init(self, context):
        self.outputs.new('RFDomainSocketType', "Meshing lines")
        self.inputs.new('NodeSocketGeometry', "Geometry Input")
        self.inputs.new('NodeSocketBool', "Vertex Group Input")

    def update(self):
        pass

    def draw_buttons(self, context, layout):
        pass

    def draw_label(self):
        return "Vertex Group to Mesh Lines"

class NODE_OT_RunSimButton(bpy.types.Operator):
    bl_idname = "node.runsim_button"
    bl_label = "run SIM"

    def execute(self, context):
        self.report({'INFO'}, "Simulation complete")
        return {'FINISHED'}
    
class NODE_OT_PreviewPECButton(bpy.types.Operator):
    bl_idname = "node.previewpec_button"
    bl_label = "Preview PEC"

    def execute(self, context):
        self.report({'INFO'}, "Simulation complete")
        return {'FINISHED'}

# Register and add to the node menu
def register():
    bpy.utils.register_class(RFDomainSocket)
    bpy.utils.register_class(MyObjectSocket)
    bpy.utils.register_class(RFSimNode)
    bpy.utils.register_class(RFSimVertexGroupToMeshLinesNode)
    bpy.utils.register_class(NODE_OT_RunSimButton)
    bpy.utils.register_class(NODE_OT_PreviewPECButton)

    def draw_node_menu(self, context):
        layout = self.layout
        layout.operator('node.add_node', text="RF Simulation").type = 'RFSimNodeType'
        layout.operator('node.add_node', text="RF Simulation").type = 'RFSimVertexGroupToMeshLinesNodeType'

    bpy.types.NODE_MT_add.append(draw_node_menu)

def unregister():
    bpy.utils.unregister_class(RFDomainSocket)
    bpy.utils.unregister_class(MyObjectSocket)
    bpy.utils.unregister_class(RFSimNode)
    bpy.utils.unregister_class(RFSimVertexGroupToMeshLinesNode)
    bpy.utils.unregister_class(NODE_OT_RunSimButton)
    bpy.utils.unregister_class(NODE_OT_PreviewPECButton)
    bpy.types.NODE_MT_add.remove(draw_node_menu)

if __name__ == "__main__":
    register()

Exploring python's API documentation can be difficult, turns out GPT-4o is quite good at that (GPT3.5 hallucinates a lot of stuff, especially the newest API accessible functions)

I haven't been able to add multi-input sockets to custom nodes though, nor attribute fields. I am not sure this is possible yet

[samerps]

Many thanks for sharing that python script, yes it might not be possible yet to add attribute fields from what I've read. I will be experimenting with these custom nodes to see if I can use them for my upcoming addon for integrating FastHenry. I'm still learning Python so it will take me some time.

For the surface current dump, I've attached my project files that I created for this animation https://www.youtube.com/watch?v=KkG2bQndp0o
In the blender file I have two current dump planes that are placed directly on the PCB copper, and in the Matlab/Octave script I have this
CSX = AddDump(CSX,'Ct_dump1', 'DumpType', 3);

PCB_ant_layoutsim.zip

[Juleinn]

Thank you very much ! :)

[Juleinn]

Some updates !

Turns out you can't run custom code on nodes dependency graph update, which would make a node-based workflow in blender for OpenEMS quite hacky and brittle.

As for displaying field dumps as volumetric data I was able to achieve some good results with a script, but integrating that into blender would probably be difficult, because of conflicting required python versions and dependencies. At this point it might be the better choice to have blender act as a network based-client, and defer all processing OpenEMS, OpenVDB, vtk, ... to a docker image as to ensure a reproducible behavior and make it easy to install.

Anyway here is the results for volumetric E/H-field rendering !

Note: what you see is not truly wavefronts, but rather an "iso - E-field", meaning it does not take into account the fact that local E-field magnitude reduces as you go further from the antenna, it just sharply displays all the points that are within a magnitude range. Displaying the whole E-field just makes a fog that fills the whole dump-box and you can barely see anything.

rotH-compress.mp4

out-compresss.mp4

and here is the script for converting vtr files to vdb :
It does an interpolation because the vtr files can have a variable-size grid, causing spaces to appear in the volumetric render. So this creates a new evenly spaced-grid of cellsize comparable to the original grid smallest cell-size (make it bigger for faster processing). Log scaling for better rendering, otherwise only the few milimeters above the antenna are visible.

import vtk
from vtk.util.numpy_support import vtk_to_numpy
import numpy as np
import scipy
import pyopenvdb as vdb
import sys
from scipy.interpolate import RegularGridInterpolator
import glob
import threading

def vtr_to_vdb(vtr_file, vdb_file):
    reader = vtk.vtkXMLRectilinearGridReader()
    reader.SetFileName(vtr_file)
    reader.Update()

    output = reader.GetOutput()

    point_data = output.GetPointData()
    a0 = point_data.GetArray(0)
    points = vtk.vtkPoints()
    output.GetPoints(points)

    point_data = vtk_to_numpy(a0)
    point_data_magnitude = np.linalg.norm(point_data, axis=1)
    points = vtk_to_numpy(points.GetData())


    #print(points)

    points_x = points[:,0]
    points_y = points[:,1]
    points_z = points[:,2]

    dx  = np.diff(points_x)
    dy  = np.diff(points_y)
    dz  = np.diff(points_z)

    dx = dx[dx != 0]
    dy = dy[dy != 0]
    dz = dz[dz != 0]

    minx = np.min(np.abs(dx))
    miny = np.min(np.abs(dy))
    minz = np.min(np.abs(dz))
    cell_size = min(minx, miny, minz)
    #print(cell_size)

# artificially force cell size up
    cell_size *= 4

# this might be off by one or two
# might also need to have a per dimension cell size
    grid_dim_x = int((np.max(points_x) - np.min(points_x)) / cell_size)
    grid_dim_y = int((np.max(points_y) - np.min(points_y)) / cell_size)
    grid_dim_z = int((np.max(points_z) - np.min(points_z)) / cell_size)

    grid_x = np.linspace(np.min(points_x), np.max(points_x), grid_dim_x)
    grid_y = np.linspace(np.min(points_y), np.max(points_y), grid_dim_y)
    grid_z = np.linspace(np.min(points_z), np.max(points_z), grid_dim_z)
    #print(f"{np.min(points_x)=}")
    #print(f"{np.min(points_y)=}")
    #print(f"{np.min(points_z)=}")
    #print(f"{np.max(points_x)=}")
    #print(f"{np.max(points_y)=}")
    #print(f"{np.max(points_z)=}")

    #print(f"target {grid_x.shape=}")
    #print(f"target {grid_y.shape=}")
    #print(f"target {grid_z.shape=}")


# Create an interpolator
    ux = np.unique(points_x)
    uy = np.unique(points_y)
    uz = np.unique(points_z)
    #print(f'input data shape{len(ux)}, {len(uy)}, {len(uz)}')
# this is the correct input data
    point_data_magnitude = point_data_magnitude.reshape(len(uz), len(uy), len(ux)).T
    interpolator = RegularGridInterpolator((ux, uy, uz), point_data_magnitude)
    #print('created interpolator')

# Generate the regular grid for interpolation
    regular_points = np.array(np.meshgrid(grid_x, grid_y, grid_z, indexing='ij')).T
    #print(f"before reshape {regular_points.shape=}")
    regular_points = regular_points.reshape(-1, 3)
    #print('Created regular points matrix')
    #print(f"after reshape {regular_points.shape=}")
#for i in range(300):
#    print(f"{regular_points[i]=}")

# Interpolation
    interpolated_values = interpolator(regular_points)
    #print('\ninterpolation complete')
    #print(f"{interpolated_values.shape=}")

# Reshape result
    interpolated_volume = interpolated_values.reshape(len(grid_z), len(grid_y), len(grid_x)).T
    # log view of the output because of massive scale between min and max values
    do_log = 1
    if do_log:
        interpolated_volume = np.log(interpolated_volume, 
                                    out=np.zeros_like(interpolated_volume), 
                                    where=interpolated_volume!=0)
        # need to add the minimimum to the log'd values because otherwise we get negative attributes
        # from values <1 pre-log
        # need the minimum accross all frames
        min_value = 15
        print(f"{min_value=}")
        interpolated_volume[interpolated_volume != 0] += min_value
        # remove negative values all together
        interpolated_volume[interpolated_volume < 0] = 0 
        # reduce the total output range
        interpolated_volume *= .01 

    #print(f"{interpolated_volume.shape=}")
    #print('complete')
    #print(interpolated_volume.shape)

### --------------------
### VDB stuff now
### --------------------

    vdb_grid = vdb.FloatGrid()
    vdb_grid.background = 0.0
    vdb_grid.name = 'temperature'

    accessor = vdb_grid.getAccessor()

# compute the scale factor 
    scale_factor = 1 / cell_size
    #print(f"scale_factor = {scale_factor} ( scale down by a factor of {cell_size})")

    for index_x, x in enumerate(grid_x):
        for index_y, y in enumerate(grid_y):
            for index_z, z in enumerate(grid_z):
                ix = int(x * scale_factor)
                iy = int(y * scale_factor)
                iz = int(z * scale_factor)
                accessor.setValueOn((ix, iy, iz),
                                    interpolated_volume[index_x][index_y][index_z])

#print(point_data_magnitude.shape)
#
#print(point_data_magnitude)
#point_data_magnitude = point_data_magnitude.reshape(len(uz), len(uy), len(ux)).T
#print(point_data_magnitude)
#
#for index_x, x in enumerate(ux):
#    for indey_y, y in enumerate(uy):
#        for indez_z, z in enumerate(uz):
#            ix = int(x * scale_factor * 2)
#            iy = int(y * scale_factor * 2)
#            iz = int(z * scale_factor * 2)
#            accessor.setValueOn((ix, iy, iz),
#                                point_data_magnitude[index_x][indey_y][indez_z])

    print(f"write {vdb_file}")
    vdb.write(vdb_file, grids=[vdb_grid])

vtr_to_vdb(sys.argv[1], sys.argv[2])

#input_glob = "E_field*vtr"
#thread_count = 18
#
#def thread_func(file_split):
#    for index, file in file_split:
#        vtr_to_vdb(file, f'{index}.vdb')
#    print('Thread complete')
#
## Inpput files have discontinuous indices, and blender tries to index 
## output vdb sequence frame on index, causing mostly empty frames
## associate a global index to be used in each thread for output file naming
#files_splits = np.array_split(np.array(list(enumerate(sorted(glob.glob("E_field*vtr"))))), thread_count)
#for file_split in files_splits:
#
#    thread = threading.Thread(target=thread_func, args=(file_split,))
#    thread.start()
#

and with this script you can run this processing on all your cores

import subprocess
import numpy as np
import glob
import threading

input_glob = "E_field*vtr"
thread_count = 18

def thread_func(file_split):
    for index, file in file_split:
        command = ['python3.10', 'vtk_note.py', file, f"{index}.vdb"]
        subprocess.run(command)
    print('Thread complete')

# Inpput files have discontinuous indices, and blender tries to index 
# output vdb sequence frame on index, causing mostly empty frames
# associate a global index to be used in each thread for output file naming
files_splits = np.array_split(np.array(list(enumerate(sorted(glob.glob("E_field*vtr"))))), thread_count)
for file_split in files_splits:

    thread = threading.Thread(target=thread_func, args=(file_split,))
    thread.start()

[samerps]

awesome, thanks for sharing. I will be experimenting with this soon

[SergeiDr]

Hi,

I saw the demo/tutorial on Youtube about your IntuitionRF plugin and would like to try it. I am not very familiar with Blender and know just little about Python. Nevertheless I tried to install your plugin and got the following error:

Traceback (most recent call last):
File "C:\Program Files\Blender Foundation\Blender 4.1\4.1\scripts\modules\addon_utils.py", line 376, in enable
mod = importlib.import_module(module_name)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Program Files\Blender Foundation\Blender 4.1\4.1\python\Lib\importlib_init_.py", line 126, in import_module
return _bootstrap._gcd_import(name[level:], package, level)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "", line 1204, in _gcd_import
File "", line 1176, in _find_and_load
File "", line 1147, in _find_and_load_unlocked
File "", line 690, in _load_unlocked
File "", line 940, in exec_module
File "", line 241, in call_with_frames_removed
File "C:\Users\MERA\AppData\Roaming\Blender Foundation\Blender\4.1\scripts\addons\IntuitionRF-master_init.py", line 13, in
from CSXCAD import CSXCAD
ModuleNotFoundError: No module named 'CSXCAD'

I have Windows 10, Blender 4.1.1 installed. OpenEMS is installed as a Ready-Made Package. I ran some examples included in openEMS and they seems to be working ok. Please help me to install your plugin correctly.

[Juleinn]

Hi,
Thank you for your interest !

There is indeed a bug, it is fixed in one of the dev branches but hasn't made it to the main one yet.
I have opened an issue, we can further discuss it there : Juleinn/IntuitionRF#1 (comment)

[mfaccin]

this is awesome guys! would it be possible to create the radiation pattern for those shown on Blender?

[Juleinn]

Yes details on that in the tutorial (2nd example)

[samerps]

@Juleinn , are you on social media, Twitter and LinkedIn?

[Juleinn]

I am on LinkedIn : linkedin.com/in/anton-claes/
I am not on Twitter but I publish blender renders on Instagram (including antennas sims) @_juleinn

[MuckYu]

Does it also take the environment materials for consideration? For example if I make a partial enclosure/wall out of metal etc.?

And can I potentially calculate/view the range/distance a signal could have?

[Juleinn]

Currently the blender plugin only supports material defined by epsilon/kappa values, and perfect electrical conductors (PEC = metals). I think OpenEMS supports anisotropic materials but it is not supported in the plugin, and I'm not sure how to do that with the python interface. Regarding walls, the simulation is defined in the plugin as having all boundaries perfectly matched layers (PML) meaning no reflections at the boundaries of the simulation. OpenEMS also supports defining PEC as boundaries, typically used to make ground planes. This is not supported in the plugin yet.

As for calculating range, OpenEMS will only give you the radiation pattern and peak gain of your antenna. To compute range you will also need to know your transmit power, receiver sensitivity and model your path loss. I believe that's out of OpenEMS' scope but there is plenty of link budget calculator utilities online where you can plug in the parameters for your system.

[istra711]

@Juleinn - Have you ever used the plugin on Linux (Ubuntu)? I tried several different distributions without success. The problem with blender is, that it is impossible to run/build the plugin in an isolated virtual python environment. The plugin has some dependencies that will break the python installation that comes with the distribution, in my case Ubuntu 22.04. Can you name a distribution that works without heavy tweaking?

[Juleinn]

Yes the plugin was actually developped with ubuntu 22.04, but the bindings were indeed brittle from the start, copying the OpenEMS sys.path is basically a hack.

Its also been tested on manjaro linux. I've seen there is discussion to make an OpenEMS python package which would simplify installation significantly.

Could you open an issue on the plugin's repo with the terminal logs from blender while using the plugin ?