This page is intended to serve as a starting point for hackers
seeking to understand, or extend, the nittygritty implementation
details of the scuffem core library.
More technical details may be found in the
libscuff Implementation Notes
and Technical Reference, available as a PDF document.
The toplevel data structure in scuffem
is a C++ class named
RWGGeometry. The definition of this class is a
little too big to present in full here (you can find it in the
header file libscuff.h ), but we will point out its
most important data fields and class methods.
Geometries in scuffem are represented
by a collection of two or more contiguous
threedimensional regions
bounded by one or more twodimensional surfaces.
Material properties (permittivity and permeability) are homogeneous
(spatially constant) in each region and described by a single
scuffem material description.
Each region is assigned an integer index starting from 0.
The RWGGeometry includes the following data fields for
identifying physical regions.
class RWGGeometry
{
...
int NumRegions;
char **RegionLabels;
MatProp **RegionMPs;
...
};
Here RegionLabels[i] is a string description for the
ith region in the problem, and RegionMPs[i] is
a pointer to an instance of MatProp describing its
frequencydependent material properies. (MatProp is
a very simple class, implemented by the
libmatprop
submodule of scuffem, for
handling frequencydependent material properties.)
RWGGeometry always starts off with a single
region (region 0 ) with label Exterior
and the material properties of vacuum.
Each REGION statement in the
.scuffgeo
file then creates a new region, starting with region 1 .
(This is true unless the label specified to the REGION
keyword is Exterior, in which case the
statement just redefines the material properties of region 0 ).
Each OBJECT...ENDOBJECT section in the .scuffgeo
file also creates a single new region (for the interior of the object).
Regions in a geometry are separated from one another by surfaces.
Each surface is described by a C++ class named RWGSurface.
The RWGGeometry class maintains an internal array of
RWGSurfaces:
class RWGGeometry
{
...
int NumSurfaces;
RWGSurface **Surfaces;
...
};
Each OBJECT...ENDOBJECT or SURFACE...ENDSURFACE
section in the .scuffgeo file adds a new RWGSurface
structure to the geometry. (Note that, in scuffem
parlance, an "object" is just a special case of a surface in which the
surface is closed.)
The RWGSurface class is again slightly
too complicated to list in full here, but we will discuss its most salient
fields and methods. Among these are the RegionIndices field:
class RWGSurface
{
...
int RegionIndices[2];
...
};
These two integers are the indices of the regions on the two sides of the
surface. The first region (RegionIndex[0] ) is
the positive region for the surface; this means that the electric
and magnetic surface currents on the surface contribute to the fields
in that region with a positive sign. The second region
(RegionIndex[1] ) is the negative region; currents
on the surface contribute to the fields in that region with a negative
sign.
Another way to think of this is that the surface normal vector
n points
away from RegionIndex[1]
and
into RegionIndex[0] .
The mesh describing each surface in a geometry is structure is analyzed
into lists of vertices, triangular panels, and panel
edges. Several internal data fields in the RWGSurface
class are devoted to storing this information.
class RWGSurface
{
...
int NumVertices;
double *Vertices;
int NumPanels;
RWGPanel **Panels;
int NumEdges;
RWGEdge **Edges;
...
};
Here Vertices is an array of 3*NumVertices
double values in which the cartesian coordinates of each
vertex are stored one after another. Thus, the x, y, z
coordinates of the nv th vertex are
Vertices[3*nv+0], Vertices[3*nv+1], Vertices[3*nv+2].
Panels and Edges are arrays of pointers
to specialized data structures for storing geometric data.
The RWGPanel and RWGEdge structures
The elemental data structure in the scuffem
geometry hierarchy is RWGPanel . Each instance of this structure
describes a single triangular panel in the mesh discretization of an
RWGSurface.
RWGPanel definition 
typedef struct RWGPanel
{
int VI[3]; /* indices of vertices in Vertices array */
int EI[3]; /* indices of edges in Edges array */
double Centroid[3]; /* panel centroid */
double ZHat[3]; /* normal vector */
double Radius; /* radius of enclosing sphere */
double Area; /* panel area */
int Index; /* index of this panel within RWGSurface */
} RWGPanel;

Here the elements of the VI array are the indices of the
three panel vertices within the list of vertices for the given
RWGSurface. The Index field in RWGPanel
indicates that panel's index within the Panels array
of the parent RWGSurface. The remaining fields tabulate some useful
geometric data on the panel.
In addition to an RWGPanel structure for each triangle
in the surface mesh, we also create an RWGEdge structure
for each panel edge.
RWGEdge definition 
typedef struct RWGEdge
{
int iV1, iV2, iQP, iQM; /* indices of panel vertices (iV1<iV2) */
double Centroid[3]; /* edge centroid */
double Length; /* length of edge */
double Radius; /* radius of enclosing sphere */
int iPPanel; /* index of PPanel within RWGSurface (0..NumPanels1)*/
int iMPanel; /* index of MPanel within RWGSurface (0..NumPanels1)*/
int PIndex; /* index of this edge within PPanel (0..2)*/
int MIndex; /* index of this edge within MPanel (0..2)*/
int Index; /* index of this edge within RWGSurface (0..NumEdges1)*/
RWGEdge *Next; /* pointer to next edge in linked list */
} RWGEdge;

The iQP, iV1, iV2, iQM fields here are indices into the
list of vertices for the parent RWGSurface.
iV1 and iV2 are the actual endpoints of the edge.
iQP and iQM denote respectively the current
source and sink vertices for the RWG basis function corresponding
to this edge.
iPPanel and iMPanel are the indices (into the
Panels array of the parent RWGSurface ) of the
positive and negative panels associated with the edge.
(The positive panel is the one from which the RWG current emanates;
its vertices are iQP, iV1, iV2.
The negative panel is the one into which the RWG current is sunk;
its vertices are iQM, iV1, iV2. )
The iPPanel and iMPanel fields take values
between 0 and NumPanels1 , where
NumPanels is defined in the parent RWGSurface.
PIndex and MIndex are the indices of the edge
within the positive and negative RWGPanel s.
(The index of an edge within a panel is defined as the index
within the panel of the panel vertex opposite that edge; the
index of a vertex is its position in the VI array in the
RWGPanel structure.)
The Index field of RWGEdge is the index of
the structure within the Edges array of the parent
RWGSurface.
The remaining fields store some geometric data on the edge.
The Centroid field stores the cartesian coordinates
of the midpoint of the line segment between vertices V1
and V2 .
Here's an example of two panels (panel indices 17 and
39 ) in a surface mesh, and a single
edge (edge index 24 ) shared between them.
Panel 17 is the positive panel for the edge,
while panel 39 is the negative panel for the edge; the corresponding
direction of current flow is indicated by the arrows.
The larger red numbers near the vertices are the indices of those
vertices in the Vertices array.
The smaller magenta and cyan numbers are the indices of the vertices
within the VI arrays in the two RWGPanels .
(The other four edges of this panel pair would also have corresponding
RWGEdge structures; these are not shown in the figure.)
The RWGGeometry class contains a hierarchy of
routines for assembling the BEM matrix. If you are simply
using the code to solve scattering problems, you will
only ever need to call the toplevel routine
AssembleBEMMatrix or perhaps the
secondhighestlevel routine, AssembleBEMMatrixBlock.
However, developers of surfaceintegralequation methods
may wish to access the lowerlevel routines for
computing the interactions of individual RWG basis functions
or for computing certain integrals over triangular regions.
AssembleBEMMatrix
The toplevel matrix assembly routine is
AssembleBEMMatrix. This routine loops over
all unique pairs of RWGSurfaces in the geometry.
For each unique pair
of surfaces, the routine calls AssembleBEMMatrixBlock
to assemble the subblock of the BEM matrix corresponding to a single
pair of RWGSurfaces, then stamps this subblock
into the appropriate place in the overall BEM matrix.
For example, if a geometry contains three RWGSurfaces ,
then its overall BEM matrix has the block structure
where the M_{ij} subblock describes the
interactions of surfaces i and j.
In this case AssembleBEMMatrix proceeds by making 6 calls
to AssembleBEMMatrixBlock , one for each of the
diagonal and abovediagonal blocks. (The belowdiagonal
blocks are related to their abovediagonal counterparts by
symmetry.)
AssembleBEMMatrixBlock
AssembleBEMMatrixBlock computes the subblock of the
BEM matrix corresponding to a single pair of RWGSurface s.
For compact (nonperiodic) geometries, this amounts to making just
a single call to SurfaceSurfaceInteractions . For
periodic geometries, this involves making multiple calls to
SurfaceSurfaceInteractions in which one of the two
surfaces is displaced through various lattice vectors L to account
for the contribution of periodic images.
The contribution of the
matrix subblock computed by SurfaceSurfaceInteractions
with displacement vector L is weighted in the overall
BEM matrix by a Bloch phase factor
e^{ik•L}
where k is the Bloch wavevector, i.e. we have
where the L superscript indicates that the corresponding
matrix subblock is to be computed with one of the two surfaces
displaced through translation vector L.
SurfaceSurfaceInteractions
SurfaceSurfaceInteractions loops over all RWGEdge
structures on each of the two surfaces it is considering.
For each pair of edges, it calls EdgeEdgeInteractions
to compute the inner products of the RWG basis functions
with the G and C dyadic Green's functions
for each of the material regions through which the two surfaces
interact. (Surfaces may interact through 0, 1, or 2 material
regions.)
Then it stamps these values into their appropriate places
in the BEM matrix subblock.
For example, if we have surfaces
S_{α} and S_{β} that
interact through a single dielectric medium, the structure of
the corresponding matrix subblock is
where e.g. b_{αm} is the
mth basis function on surface S_{α}.
The kernels
here are scalar multiples of the G and C
dyadics:
where k, Z are the wavenumber and (absolute)
wave impedance of the dielectric medium at the frequency
in question.
The matrix structure above is for the case of two surfaces
interacting through a single material region (for example,
S_{α} and S_{β}
might be the outer surfaces of two compact objects embedded
in vacuum or in a homogeneous medium, in which case the
surfaces interact only throught that medium).
If the surfaces interact through two material regions
(for example, if we have
S_{α}=S_{β}
and we are computing the selfinteraction of the outer
surface of a dielectric object embedded in vacuum) then
each matrix entry is actually a sum of two
inner products, one involving the Γ kernels for
the interior medium and one involving the kernels for
the exterior. (If the surfaces are PEC, then the dimension
of the matrix is halved with only the EE terms
retained.)
EdgeEdgeInteractions
EdgeEdgeInteractions considers a pair of RWG basis
functions and computes the inner products of these basis functions
with the G and C dyadic Green's functions for a
single material medium. Because each basis function is supported
on two triangles, the full inner products involve sums of four
trianglepair contributions:
Here l_{m},l_{n} are the lengths of the
interior edges to which the RWG basis functions are associated, and
each of the four terms in the sum is a fourdimensional integral
over a single pair of panels, computed by
PanelPanelInteractions.
PanelPanelInteractions
PanelPanelInteractions is the lowestlevel routine in the
scuffem BEM matrix assembly hierarchy.
This routine computes the individual terms in equation (1) above,
i.e. the contributions of a single pair of panels
to the inner products of two RWG basis functions with the G and
C dyadic Green's functions:
Note that the RWG basisfunction prefactor l/2A is broken
up into two factors between equations (1) and (2): the 1/2A
part is included in the panelpanel integrals (2), while
the l part is included when summing the four panelpanel
contributions in equation (1) to obtain the overall inner product.
If the panels are far away from each other, PanelPanelInteractions
uses loworder fourdimensional numerical cubature to compute the full
interactions. Otherwise, PanelPanelInteractions uses
loworder fourdimensional numerical cubature to compute the interactions
with singularitysubtracted versions of the G and C
kernels, then adds the contributions of the singular terms after
looking them up in an internallystored cache.
Here's a worked example of a matrixelement computation in
a scuffem run.
The geometry and the labeling of panels, edges, and vertices
We'll consider a
scattering geometry consisting of a single cube of dielectric
material, with side length L=1 μm, discretized into
triangles with minor side length L/10, yielding a total
of 1200 triangles, 1800 interior edges, and 3600 RWG basis
functions (1800 each for electric and magnetic currents)
for a dielectric geometry. The
gmsh geometry file for this example is
Cube_N.geo,
the gmsh mesh file produced by running
gmsh 2 Cube_N.geo is
Cube_10.msh,
and a scuffem geometry file describing
a geometry consisting of this discretized cube with interior dielectric
permittivity ε=4 is
DielectricCube.scuffgeo.
(Click the image below for higher resolution.)
For lowlevel work it is convenient to have data on the internal
indices that scuffem assigns to
panels, vertices, and edges when it reads in a geometry file.
We get this by running scuffanalyze on the mesh file
in question with the WriteGMSHLabels options:
% scuffanalyze mesh Cube_10.msh WriteGMSHFiles WriteGMSHLabels
This produces a file named Cube_10.pp , which we open in
gmsh to produce a graphical depiction of
the labeling of panels, edges, and vertices. Zooming in on the region
near the origin, we will focus on the two panels lying closest to the
origin in the xy plane; we see that
scuffem has assigned these two panels
panel indices 0 and 1 , respectively, while
the edge they share has been assigned interioredge index 0 .
RWG basis function b_{0}
Here's a schematic depiction of the panels that comprise the
basis function b_{0} associated with interior edge
0 :
In this diagram, the edge length is L=0.1, O denotes
the origin of coordinates, and the vertices marked
Q^{±} are the
source and sink vertices for the current distribution.
The RWG basis function b_{0} associated with
interior edge 0 is
where the RWG basis function edge length is l=√2L
and the panel areas are
A_{0}=A_{1}=L^{2}/2.
Note that
P_{0} and
P_{1} are respectively the positive and negative
panels associated with basis function b_{0}.
Panelpanel interactions
Here's a C++ code snippet that computes the
quantities
G_{00}^{++}
and
G_{00}^{+} in equation (1a)
above, i.e. the contributions of the positivepositive
and positivenegative panel pairs to the inner product of
b_{0} with itself through the G kernel,
with the wavenumber set to k=1.0;
Just to be totally explicit, the numbers that are being
computed here are
// read in geometry from .scuffgeo file
RWGGeometry *G = new RWGGeometry("DielectricCube.scuffgeo");
RWGSurface *S = G>Surfaces[0];
// initialize an argument structure for GetPanelPanelInteractions
GetPPIArgStruct MyPPIArgs, *PPIArgs=&MyPPIArgs;
InitGetPPIArgs(PPIArgs);
PPIArgs>Sa = PPIArgs>Sb = S;
PPIArgs>k = 1.0;
// fill in arguments to compute the positivepositive panel pair
PPIArgs>npa = 0;
PPIArgs>iQa = 1;
PPIArgs>npb = 0;
PPIArgs>iQb = 1;
GetPanelPanelInteractions(PPIArgs);
printf("G++ = %e + %ei\n",real(PPIArgs>H[0]),imag(PPIArgs>H[0]));
// fill in arguments to compute the positivenegative panel pair
PPIArgs>npa = 0;
PPIArgs>iQa = 1;
PPIArgs>npb = 1;
PPIArgs>iQb = 2;
GetPanelPanelInteractions(PPIArgs);
printf("G++ = %e + %ei\n",real(PPIArgs>H[0]),imag(PPIArgs>H[0]));
This code produces the following output:
G++ = 3.189105e+00 + 7.950381e02i
G+ = 1.537561e+00 + 7.956272e02i
Notice that, to compute a panelpanel interaction, you allocate
and initialize an instance of a data structure called
GetPPIArgs . This structure contains a large
number of fields which in many cases can be set to default
values; to ensure that these fields are properly initialized,
always call InitGetPPIArgs() on an newlyallocated
instance of GetPPIArgs .
Then, you fill in the appropriate fields of this structure
to specify the two panels over which you want to integrate.
Specifically, the fields
Sa , npa
specify the first panel (the npa th RWGPanel
in the Panels array for surface Sa ),
while the field iQa (an integer in the range 0..2 )
identifies the index of the Q vertex (RWG current source/sink
vertex) within the three vertices of the panel.
The fields
Sb , npb , and iQb
similarly identify the second panel and source/sink vertex.
Initialize the k field in the
GetPPIArgs structure to the wavenumber parameter
in the Helmholtz kernel. (k may be complex or purely imaginary.)
Edgepanel interactions
Here's a C++ code snippet that computes the
full basisfunction inner product
< b_{0}  G  b_{0} >
in equation (1a) above.
GetEEIArgStruct MyEEIArgs, *EEIArgs=&MyEEIArgs;
InitGetEEIArgs(EEIArgs);
EEIArgs>Sa = EEIArgs>Sb = S;
EEIArgs>nea = EEIArgs>neb = 0;
EEIArgs>k = 1.0;
GetEdgeEdgeInteractions(EEIArgs);
printf("<bGb> = %e + %ei\n",real(EEIArgs>GC[0]),imag(EEIArgs>GC[0]));
Note that, similar to the case of GetPanelPanelInteractions
discussed above, the call to GetEdgeEdgeInteractions takes as
argument a pointer to a struct of type GetEEIArgStruct .
As before, you should always call InitGetEEIArgs
to initialize this structure, then set whichever fields you
need to specify.
In this case, the only fields we need to set are
Sa,Sb (indices of the RWGSurface ),
nea,neb
(indices of the interior edges in the Edges array
corresponding to the RWG basis functions), and k
(wavenumber parameter in the Helmholtz kernel).
The result of this code snippet is
<bGb> = 6.606176e02 + 2.356554e06i
Using equation (1a) above, we can understand this result in conjunction
with the results printed out above for G++ and
G+ . For this particular basis function, only two
of the four panelpanel interactions on the RHS of equation (1a)
have distinct values (because the two panels that comprise the
basis function have the same shapes and areas, so we have
G+ = G+ and G = G++ ), and the
edgelength prefactors l_{m}, l^{n}
both have value 0.1•√2. Thus for this case
equation (1a) reads
< b_{0}  G  b_{0} >
=2•0.02•(G++  G+ )
and, indeed, plugging in the numbers printed out above, we find
{6.6e02, 2.4e06i}
= 2•0.02•( {3.19,7.95e2i}  {1.54,7.96e2i} ).
