# Data Structures and Class Methods in scuff-em

This page is intended to serve as a starting point for hackers seeking to understand, or extend, the nitty-gritty implementation details of the scuff-em core library.

More technical details may be found in the libscuff Implementation Notes and Technical Reference, available as a PDF document.

- Geometries, Regions, Surfaces
- Panels, Edges, Vertices
- Assembling the BEM Matrix
- An Explicit Low-Level Example

** Data Structures and Class Methods in scuff-em**

## 1. Geometries, Regions, Surfaces

The top-level data structure in libscuff
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 scuff-em are represented
a collection of two or more contiguous
three-dimensional *regions*
bounded by one or more two-dimensional *surfaces.*
Material properties (permittivity and permeability) are homogeneous
(spatially constant) in each region and described by a single
scuff-em 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
*i*th region in the problem, and `RegionMPs[i]`

is
a pointer to an instance of `MatProp`

describing its
frequency-dependent material properies. (`MatProp`

is
a very simple class, implemented by the
libmatprop submodule of scuff-em, for
handling frequency-dependent 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 scuff-em
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]`

.

## 2. Panels, Edges, Vertices

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 scuff-em
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..NumPanels-1)*/ int iMPanel; /* index of MPanel within RWGSurface (0..NumPanels-1)*/ 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..NumEdges-1)*/ 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 `NumPanels-1`

, 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.)

## 3. Assembling the BEM Matrix

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 top-level routine
`AssembleBEMMatrix`

or perhaps the
second-highest-level routine, `AssembleBEMMatrixBlock.`

However, developers of surface-integral-equation methods
may wish to access the lower-level routines for
computing the interactions of individual RWG basis functions
or for computing certain integrals over triangular regions.

`AssembleBEMMatrix`

The top-level 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 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 above-diagonal blocks. (The below-diagonal
blocks are related to their above-diagonal counterparts by
symmetry.)

`AssembleBEMMatrixBlock`

`AssembleBEMMatrixBlock`

computes the subblock of the
BEM matrix corresponding to a single pair of `RWGSurface`

s.
For compact (non-periodic) 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
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
*m*th 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 self-interaction 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
triangle-pair 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 four-dimensional integral
over a single pair of panels, computed by

`PanelPanelInteractions.`

`PanelPanelInteractions`

`PanelPanelInteractions`

is the lowest-level routine in the
scuff-em 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 basis-function prefactor *l/2A* is broken
up into two factors between equations (1) and (2): the *1/2A*
part is included in the panel-panel integrals (2), while
the *l* part is included when summing the four panel-panel
contributions in equation (1) to obtain the overall inner product.

If the panels are far away from each other, `PanelPanelInteractions`

uses low-order four-dimensional numerical cubature to compute the full
interactions. Otherwise, `PanelPanelInteractions`

uses
low-order four-dimensional numerical cubature to compute the interactions
with *singularity-subtracted versions* of the **G** and **C**
kernels, then adds the contributions of the singular terms after
looking them up in an internally-stored cache.

## 4. An Explicit Low-Level Example

Here's a worked example of a matrix-element computation in a scuff-em 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 scuff-em geometry file describing
a geometry consisting of this discretized cube with interior dielectric
permittivity ε=4 is
DielectricCube.scuffgeo.

For low-level work it is convenient to have data on the internal
indices that scuff-em assigns to
panels, vertices, and edges when it reads in a geometry file.
We get this by running `scuff-analyze`

on the mesh file
in question with the `--WriteGMSHLabels`

options:

% scuff-analyze --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
scuff-em has assigned these two panels
panel indices `0`

and `1`

, respectively, while
the edge they share has been assigned interior-edge index `0`

.

## RWG basis function **b**_{0}

_{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*=√2*L*
and the panel areas are
*A _{0}=A_{1}*=

*L*

^{2}/2. Note that

*P*and

_{0}*P*are respectively the positive and negative panels associated with basis function

_{1}**b**

_{0}.

### Panel-panel interactions

Here's a `C++`

code snippet that computes the
quantities
*G _{00}^{++}*
and

*G*in equation (1a) above, i.e. the contributions of the positive-positive and positive-negative panel pairs to the inner product of

_{00}^{+-}**b**with itself through the

_{0}**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 positive-positive 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 positive-negative 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.950381e-02i G+- = -1.537561e+00 + -7.956272e-02i

Notice that, to compute a panel-panel 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 newly-allocated
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.)

### Edge-panel interactions

Here's a `C++`

code snippet that computes the
full basis-function 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("<b|G|b> = %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

<b|G|b> = -6.606176e-02 + 2.356554e-06i

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 panel-panel 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

edge-length 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.6e-02, 2.4e-06i} = 2•0.02•( {-3.19,-7.95e-2i} - {-1.54,7.96e-2i} ).