python: update ports and nf2ff handling, update Tutorials

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

python/Tutorials/Helical_Antenna.py

@@ -116,6 +116,8 @@ nf2ff = FDTD.CreateNF2FFBox(opt_resolution=[lambda0/15]*3)
 
 if 0:  # debugging only
     CSX_file = os.path.join(Sim_Path, 'helix.xml')
+    if not os.path.exists(Sim_Path):
+        os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
     os.system(r'AppCSXCAD "{}"'.format(CSX_file))
 
@@ -153,25 +155,25 @@ theta = arange(0.,180.,1.)
 phi = arange(-180,180,2)
 disp( 'calculating the 3D far field...' )
 
-nf2ff.CalcNF2FF(Sim_Path, f0, theta, phi, read_cached=True, verbose=True )
+nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f0, theta, phi, read_cached=True, verbose=True )
 
 #
-Dmax_dB = 10*log10(nf2ff.Dmax[0])
+Dmax_dB = 10*log10(nf2ff_res.Dmax[0])
-E_norm = 20.0*log10(nf2ff.E_norm[0]/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.Dmax[0])
+E_norm = 20.0*log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
 
 theta_HPBW = theta[ np.where(squeeze(E_norm[:,phi==0])<Dmax_dB-3)[0][0] ]
 #
 # display power and directivity
-print('radiated power: Prad = {} W'.format(nf2ff.Prad[0]))
+print('radiated power: Prad = {} W'.format(nf2ff_res.Prad[0]))
 print('directivity: Dmax = {} dBi'.format(Dmax_dB))
-print('efficiency: nu_rad = {} %'.format(100*nf2ff.Prad[0]/interp(f0, freq, port.P_acc)))
+print('efficiency: nu_rad = {} %'.format(100*nf2ff_res.Prad[0]/interp(f0, freq, port.P_acc)))
 print('theta_HPBW = {} °'.format(theta_HPBW))
 
 
 ##
-E_norm = 20.0*log10(nf2ff.E_norm[0]/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.Dmax[0])
+E_norm = 20.0*log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
-E_CPRH = 20.0*log10(np.abs(nf2ff.E_cprh[0])/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.Dmax[0])
+E_CPRH = 20.0*log10(np.abs(nf2ff_res.E_cprh[0])/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
-E_CPLH = 20.0*log10(np.abs(nf2ff.E_cplh[0])/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.Dmax[0])
+E_CPLH = 20.0*log10(np.abs(nf2ff_res.E_cplh[0])/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
 
 ##
 figure()
@@ -181,7 +183,7 @@ plot(theta, E_CPLH[:,phi==0],'r-.', linewidth=2, label='$|E_{CPLH}|$')
 grid()
 xlabel('theta (deg)')
 ylabel('directivity (dBi)')
-title('Frequency: {} GHz'.format(nf2ff.freq[0]/1e9))
+title('Frequency: {} GHz'.format(nf2ff_res.freq[0]/1e9))
 legend()
 
 ### dump to vtk

python/Tutorials/MSL_NotchFilter.py

@@ -85,13 +85,16 @@ port[1] = FDTD.AddMSLPort( 2, pec, portstart, portstop, 'x', 'z', MeasPlaneShift
 ## Filter-stub
 start = [-MSL_width/2,  MSL_width/2, substrate_thickness]
 stop  = [ MSL_width/2,  MSL_width/2+stub_length, substrate_thickness]
-CSX.AddBox( pec, start, stop, priority=10 )
+box = CSX.AddBox( pec, start, stop, priority=10 )
 
 if 0:  # debugging only
     CSX_file = os.path.join(Sim_Path, 'notch.xml')
+    if not os.path.exists(Sim_Path):
+        os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
     os.system(r'AppCSXCAD "{}"'.format(CSX_file))
 
+
 if not post_proc_only:
     FDTD.Run(Sim_Path, verbose=3, cleanup=True)
 

python/Tutorials/Rect_Waveguide.py

@@ -86,6 +86,8 @@ CSX.AddBox(Et, start, stop);
 ## Write openEMS compatoble xml-file ######################################
 if 0:  # debugging only
     CSX_file = os.path.join(Sim_Path, 'rect_wg.xml')
+    if not os.path.exists(Sim_Path):
+        os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
     os.system(r'AppCSXCAD "{}"'.format(CSX_file))
 

python/Tutorials/Simple_Patch_Antenna.py

@@ -48,23 +48,20 @@ mesh_res = C0/(f0+fc)/1e-3/20
 
 CSX = CSXCAD.ContinuousStructure()
 FDTD.SetCSX(CSX)
-
+mesh = CSX.GetGrid()
+mesh.SetDeltaUnit(1e-3)
 
 #initialize the mesh with the "air-box" dimensions
-mesh = {}
+mesh.AddLine('x', [-SimBox[0]/2, feed_pos, SimBox[0]/2])
-mesh['x'] = [-SimBox[0]/2, feed_pos, SimBox[0]/2]
+mesh.AddLine('y', [-SimBox[1]/2, SimBox[1]/2]          )
-mesh['y'] = [-SimBox[1]/2, SimBox[1]/2]
+mesh.AddLine('z', [-SimBox[2]/3, SimBox[2]*2/3]        )
-mesh['z'] = [-SimBox[2]/3, SimBox[2]*2/3]
 
 ## create patch
 patch = CSX.AddMetal( 'patch' ) # create a perfect electric conductor (PEC)
 start = [-patch_width/2, -patch_length/2, substrate_thickness]
 stop  = [ patch_width/2 , patch_length/2, substrate_thickness]
 pb=CSX.AddBox(patch, priority=10, start=start, stop=stop) # add a box-primitive to the metal property 'patch'
-
+pb.AddEdges2Grid('xy', metal_edge_res=mesh_res/2)
-edge_mesh = np.array([-1/3.0, 2/3.0]) * mesh_res/2
-mesh['x'] = r_[mesh['x'], start[0]-edge_mesh, stop[0]+edge_mesh]
-mesh['y'] = r_[mesh['y'], start[1]-edge_mesh, stop[1]+edge_mesh]
 
 ## create substrate
 substrate = CSX.AddMaterial( 'substrate', epsilon=substrate_epsR, kappa=substrate_kappa)
@@ -73,31 +70,27 @@ stop  = [ substrate_width/2,  substrate_length/2, substrate_thickness]
 sb=CSX.AddBox( substrate, priority=0, start=start, stop=stop )
 
 # add extra cells to discretize the substrate thickness
-mesh['z'] = r_[mesh['z'], linspace(0,substrate_thickness,substrate_cells+1)]
+mesh.AddLine('z', linspace(0,substrate_thickness,substrate_cells+1))
 
 ## create ground (same size as substrate)
 gnd = CSX.AddMetal( 'gnd' ) # create a perfect electric conductor (PEC)
 start[2]=0
 stop[2] =0
-gb=CSX.AddBox(gnd, start, stop, priority=10)
+gb=CSX.AddBox(gnd, start, stop, priority=10, edges2grid='all')
-
-mesh['x'] = r_[mesh['x'], start[0], stop[0]]
-mesh['y'] = r_[mesh['y'], start[1], stop[1]]
 
 ## apply the excitation & resist as a current source
 start = [feed_pos, 0, 0]
 stop  = [feed_pos, 0, substrate_thickness]
-port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'z', 1.0, priority=5)
+port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'z', 1.0, priority=5, edges2grid='all')
 
-mesh['x'] = r_[mesh['x'], start[0]]
+mesh.SmoothMeshLines('all', mesh_res, 1.4)
-mesh['y'] = r_[mesh['y'], start[1]]
-
-CSX.DefineGrid(mesh, unit=1e-3, smooth_mesh_res=mesh_res)
 
 nf2ff = FDTD.CreateNF2FFBox()
 
 if 0:  # debugging only
     CSX_file = os.path.join(Sim_Path, 'simp_patch.xml')
+    if not os.path.exists(Sim_Path):
+        os.mkdir(Sim_Path)
     CSX.Write2XML(CSX_file)
     os.system(r'AppCSXCAD "{}"'.format(CSX_file))
 
@@ -121,10 +114,10 @@ else:
     f_res = f[idx[0]]
     theta = np.arange(-180.0, 180.0, 2.0)
     phi   = [0., 90.]
-    nf2ff.CalcNF2FF(Sim_Path, f_res, theta, phi, center=[0,0,1e-3], read_cached=False )
+    nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, phi, center=[0,0,1e-3])
 
     figure()
-    E_norm = 20.0*np.log10(nf2ff.E_norm[0]/np.max(nf2ff.E_norm[0])) + nf2ff.Dmax[0]
+    E_norm = 20.0*np.log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + nf2ff_res.Dmax[0]
     plot(theta, np.squeeze(E_norm[:,0]), label='xz-plane')
     plot(theta, np.squeeze(E_norm[:,1]), label='yz-plane')
     grid()

python/openEMS/nf2ff.py

@@ -69,53 +69,60 @@ class nf2ff:
                 CSX.AddBox(self.e_dump, l_start, l_stop)
                 CSX.AddBox(self.h_dump, l_start, l_stop)
 
-    def CalcNF2FF(self, sim_path, freq, theta, phi, center=[0,0,0], read_cached=True, verbose=0):
+    def CalcNF2FF(self, sim_path, freq, theta, phi, center=[0,0,0], outfile=None, read_cached=True, verbose=0):
-        if not hasattr(freq, "__iter__"):
+        if np.isscalar(freq):
             freq = [freq]
         self.freq  = freq
+        if np.isscalar(theta):
+            theta = [theta]
         self.theta = theta
+        if np.isscalar(phi):
+            phi = [phi]
         self.phi   = phi
         self.center = center
 
-        fn = os.path.join(sim_path, self.name + '.h5')
+        if outfile is None:
-        if os.path.exists(fn) and read_cached:
+            fn = os.path.join(sim_path, self.name + '.h5')
-            self.ReadNF2FF(sim_path)
+        else:
-            return
+            fn = os.path.join(sim_path,  outfile)
+        if  not read_cached or not os.path.exists(fn):
+            nfc = _nf2ff._nf2ff(self.freq, np.deg2rad(theta), np.deg2rad(phi), center, verbose=verbose)
 
-        nfc = _nf2ff._nf2ff(self.freq, np.deg2rad(theta), np.deg2rad(phi), center, verbose=verbose)
+            for ny in range(3):
+                nfc.SetMirror(self.mirror[2*ny]  , ny, self.start[ny])
+                nfc.SetMirror(self.mirror[2*ny+1], ny, self.stop[ny])
 
-        for ny in range(3):
+            for n in range(6):
-            nfc.SetMirror(self.mirror[2*ny]  , ny, self.start[ny])
+                fn_e = os.path.join(sim_path, self.e_file + '_{}.h5'.format(n))
-            nfc.SetMirror(self.mirror[2*ny+1], ny, self.stop[ny])
+                fn_h = os.path.join(sim_path, self.h_file + '_{}.h5'.format(n))
+                if os.path.exists(fn_e) and os.path.exists(fn_h):
+                    assert nfc.AnalyseFile(fn_e, fn_h)
 
-        for n in range(6):
+            nfc.Write2HDF5(fn)
-            fn_e = os.path.join(sim_path, self.e_file + '_{}.h5'.format(n))
-            fn_h = os.path.join(sim_path, self.h_file + '_{}.h5'.format(n))
-            if os.path.exists(fn_e) and os.path.exists(fn_h):
-                assert nfc.AnalyseFile(fn_e, fn_h)
 
-        nfc.Write2HDF5(fn)
+        result = nf2ff_results(fn)
+        if result.phi is not None:
+            assert utilities.Check_Array_Equal(np.rad2deg(result.theta), self.theta, 1e-4)
+            assert utilities.Check_Array_Equal(np.rad2deg(result.phi), self.phi, 1e-4)
+            assert utilities.Check_Array_Equal(result.freq, self.freq, 1e-6, relative=True)
+        return result
 
-        self.ReadNF2FF(sim_path)
+class nf2ff_results:
-
+    def __init__(self, fn):
-    def ReadNF2FF(self, sim_path):
+        self.fn  = fn
-        h5_file = h5py.File(os.path.join(sim_path, self.name + '.h5'), 'r')
+        h5_file  = h5py.File(fn, 'r')
         mesh_grp = h5_file['Mesh']
-        phi   = np.array(mesh_grp['phi'])
+        self.phi   = np.array(mesh_grp['phi'])
-        theta = np.array(mesh_grp['theta'])
+        self.theta = np.array(mesh_grp['theta'])
         self.r     = np.array(mesh_grp['r'])
 
         data  = h5_file['nf2ff']
-        freq = np.array(data.attrs['Frequency'])
+        self.freq = np.array(data.attrs['Frequency'])
 
-        if self.phi is not None:
-            assert utilities.Check_Array_Equal(np.rad2deg(theta), self.theta, 1e-4)
-            assert utilities.Check_Array_Equal(np.rad2deg(phi), self.phi, 1e-4)
-            assert utilities.Check_Array_Equal(freq, self.freq, 1e-6, relative=True)
         self.Dmax = np.array(data.attrs['Dmax'])
         self.Prad = np.array(data.attrs['Prad'])
 
-        THETA, PHI = np.meshgrid(theta, phi, indexing='ij')
+        THETA, PHI = np.meshgrid(self.theta, self.phi, indexing='ij')
         cos_phi = np.cos(PHI)
         sin_phi = np.sin(PHI)
 
@@ -125,7 +132,7 @@ class nf2ff:
         self.E_norm  = []
         self.E_cprh  = []
         self.E_cplh  = []
-        for n in range(len(freq)):
+        for n in range(len(self.freq)):
             E_theta = np.array(h5_file['/nf2ff/E_theta/FD/f{}_real'.format(n)]) + 1j*np.array(h5_file['/nf2ff/E_theta/FD/f{}_imag'.format(n)])
             E_theta = np.swapaxes(E_theta, 0, 1)
             E_phi   = np.array(h5_file['/nf2ff/E_phi/FD/f{}_real'.format(n)])   + 1j*np.array(h5_file['/nf2ff/E_phi/FD/f{}_imag'.format(n)])

python/openEMS/openEMS.pyx

@@ -11,46 +11,58 @@ cimport openEMS
 from . import ports, nf2ff
 
 cdef class openEMS:
+    """ openEMS
+
+    :param NrTS:           max. number of timesteps to simulate (e.g. default=1e9)
+    :param EndCriteria:    end criteria, e.g. 1e-5, simulations stops if energy has decayed by this value (<1e-4 is recommended, default=1e-5)
+    :param MaxTime:        max. real time in seconds to simulate
+    :param OverSampling:   nyquist oversampling of time domain dumps
+    :param CoordSystem:    choose coordinate system (0 Cartesian, 1 Cylindrical)
+    :param MultiGrid:      define a cylindrical sub-grid radius ( not implemented yet )
+    :param TimeStep:       force to use a given timestep (dangerous!)
+    :param TimeStepFactor: reduce the timestep by a given factor (>0 to <=1)
+    :param TimeStepMethod: 1 or 3 chose timestep method (1=CFL, 3=Rennigs (default))
+    :param CellConstantMaterial: set to 1 to assume a material is constant inside a cell (material probing in cell center)
+    """
     @staticmethod
     def WelcomeScreen():
         _openEMS.WelcomeScreen()
 
-    """ __cinit__
-    NrTS:           max. number of timesteps to simulate (e.g. default=1e9)
-    EndCriteria:    end criteria, e.g. 1e-5, simulations stops if energy has decayed by this value (<1e-4 is recommended, default=1e-5)
-    MaxTime:        max. real time in seconds to simulate
-    OverSampling:   nyquist oversampling of time domain dumps
-    CoordSystem:    choose coordinate system (0 Cartesian, 1 Cylindrical)
-    #MultiGrid:      define a cylindrical sub-grid radius
-    TimeStep:       force to use a given timestep (dangerous!)
-    TimeStepFactor: reduce the timestep by a given factor (>0 to <=1)
-    TimeStepMethod: 1 or 3 chose timestep method (1=CFL, 3=Rennigs (default))
-    CellConstantMaterial: set to 1 to assume a material is constant inside a cell (material probing in cell center)
-    """
     def __cinit__(self, *args, **kw):
         self.thisptr = new _openEMS()
         self.CSX = None
 
         if 'NrTS' in kw:
             self.SetNumberOfTimeSteps(kw['NrTS'])
+            del kw['NrTS']
         else:
             self.SetNumberOfTimeSteps(1e9)
         if 'EndCriteria' in kw:
             self.SetEndCriteria(kw['EndCriteria'])
+            del kw['EndCriteria']
         if 'MaxTime' in kw:
             self.SetMaxTime(kw['MaxTime'])
+            del kw['MaxTime']
         if 'OverSampling' in kw:
             self.SetOverSampling(kw['OverSampling'])
+            del kw['OverSampling']
         if 'CoordSystem' in kw:
             self.SetCoordSystem(kw['CoordSystem'])
+            del kw['CoordSystem']
         if 'TimeStep' in kw:
             self.SetTimeStep(kw['TimeStep'])
+            del kw['TimeStep']
         if 'TimeStepFactor' in kw:
             self.SetTimeStepFactor(kw['TimeStepFactor'])
+            del kw['TimeStepFactor']
         if 'TimeStepMethod' in kw:
             self.SetTimeStepMethod(kw['TimeStepMethod'])
+            del kw['TimeStepMethod']
         if 'CellConstantMaterial' in kw:
             self.SetCellConstantMaterial(kw['CellConstantMaterial'])
+            del kw['CellConstantMaterial']
+
+        assert len(kw)==0, 'Unknown keyword arguments: "{}"'.format(kw)
 
     def __dealloc__(self):
         del self.thisptr
@@ -58,41 +70,68 @@ cdef class openEMS:
             self.CSX.thisptr = NULL
 
     def SetNumberOfTimeSteps(self, val):
+        """ SetNumberOfTimeSteps(val)
+        """
         self.thisptr.SetNumberOfTimeSteps(val)
 
     def SetEndCriteria(self, val):
+        """ SetEndCriteria(val)
+        """
         self.thisptr.SetEndCriteria(val)
 
     def SetOverSampling(self, val):
+        """ SetOverSampling(val)
+        """
         self.thisptr.SetOverSampling(val)
 
     def SetCellConstantMaterial(self, val):
+        """ SetCellConstantMaterial(val)
+        """
         self.thisptr.SetCellConstantMaterial(val)
 
     def SetCoordSystem(self, val):
-        assert (val==0 or val==1)
+        """ SetCoordSystem(val)
+        """
+        assert (val==0 or val==1), 'SetCoordSystem: Invalid coordinate system'
         if val==0:
             pass
         elif val==1:
             self.SetCylinderCoords()
 
     def SetCylinderCoords(self):
+        """ SetCylinderCoords()
+        """
         self.thisptr.SetCylinderCoords(True)
 
     def SetTimeStepMethod(self, val):
+        """ SetTimeStepMethod(val)
+        """
         self.thisptr.SetTimeStepMethod(val)
+
     def SetTimeStep(self, val):
+        """ SetTimeStep(val)
+        """
         self.thisptr.SetTimeStep(val)
+
     def SetTimeStepFactor(self, val):
+        """ SetTimeStepFactor(val)
+        """
         self.thisptr.SetTimeStepFactor(val)
+
     def SetMaxTime(self, val):
+        """ SetMaxTime(val)
+        """
         self.thisptr.SetMaxTime(val)
 
     def SetGaussExcite(self, f0, fc):
+        """ SetGaussExcite(f0, fc)
+        """
         self.thisptr.SetGaussExcite(f0, fc)
 
 
     def SetBoundaryCond(self, BC):
+        """ SetBoundaryCond(BC)
+        """
         assert len(BC)==6
         for n in range(len(BC)):
             if type(BC[n])==int:
@@ -108,19 +147,27 @@ cdef class openEMS:
             raise Exception('Unknown boundary condition')
 
     def AddLumpedPort(self, port_nr, R, start, stop, p_dir, excite=0, **kw):
-        assert self.CSX is not None
+        """ AddLumpedPort(port_nr, R, start, stop, p_dir, excite=0, **kw)
+        """
+        assert self.CSX is not None, 'AddLumpedPort: CSX is not set!'
         return ports.LumpedPort(self.CSX, port_nr, R, start, stop, p_dir, excite, **kw)
 
     def AddRectWaveGuidePort(self, port_nr, start, stop, p_dir, a, b, mode_name, excite=0, **kw):
-        assert self.CSX is not None
+        """ AddRectWaveGuidePort(port_nr, start, stop, p_dir, a, b, mode_name, excite=0, **kw)
+        """
+        assert self.CSX is not None, 'AddRectWaveGuidePort: CSX is not set!'
         return ports.RectWGPort(self.CSX, port_nr, start, stop, p_dir, a, b, mode_name, excite, **kw)
 
     def AddMSLPort(self, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw):
-        assert self.CSX is not None
+        """ AddMSLPort(port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw)
+        """
+        assert self.CSX is not None, 'AddMSLPort: CSX is not set!'
         return ports.MSLPort(self.CSX, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite, **kw)
 
     def CreateNF2FFBox(self, name='nf2ff', start=None, stop=None, **kw):
-        assert self.CSX is not None
+        """ CreateNF2FFBox(name='nf2ff', start=None, stop=None, **kw)
+        """
+        assert self.CSX is not None, 'CreateNF2FFBox: CSX is not set!'
         directions = [True]*6
         mirror     = [0]*6
         BC_size = [0]*6
@@ -152,10 +199,14 @@ cdef class openEMS:
         return nf2ff.nf2ff(self.CSX, name, start, stop, directions=directions, mirror=mirror, **kw)
 
     def SetCSX(self, ContinuousStructure CSX):
+        """ SetCSX(CSX)
+        """
         self.CSX = CSX
         self.thisptr.SetCSX(CSX.thisptr)
 
     def Run(self, sim_path, cleanup=False, setup_only=False, verbose=None):
+        """ Run(sim_path, cleanup=False, setup_only=False, verbose=None)
+        """
         if cleanup and os.path.exists(sim_path):
             shutil.rmtree(sim_path)
             os.mkdir(sim_path)

python/openEMS/ports.py

@@ -116,14 +116,14 @@ class LumpedPort(Port):
         self.exc_ny  = CheckNyDir(exc_dir)
 
         self.direction = np.sign(self.stop[self.exc_ny]-self.start[self.exc_ny])
-        assert self.start[self.exc_ny]!=self.stop[self.exc_ny]
+        assert self.start[self.exc_ny]!=self.stop[self.exc_ny], 'LumpedPort: start and stop may not be identical in excitation direction'
 
         if self.R > 0:
             lumped_R = CSX.AddLumpedElement(self.lbl_temp.format('resist'), ny=self.exc_ny, caps=True, R=self.R)
         elif self.R==0:
             lumped_R = CSX.AddMetal(self.lbl_temp.format('resist'))
 
-        CSX.AddBox(lumped_R, self.start, self.stop, priority=self.priority)
+        CSX.AddBox(lumped_R, self.start, self.stop, priority=self.priority, edges2grid=kw.get('edges2grid', None))
 
         if excite!=0:
             exc_vec = np.zeros(3)