M. T. Homer Reid MIT Home Page
Physics Problems Research teaching My Music About Me Miscellany


libscuff API Documentation:
Ancillary Operations on RWGGeometry objects

This page documents the portion of the libscuff API that handles ancillary RWGGeometry operations that lie outside the main flow for scattering problems.

Ancillary RWGGeometry Operations
1. Visualization
2. Geometrical Transformations
3. Setting material properties
4. Specialized routines for computing field data

5. Routines for evaluating numerical cubatures over panels

1. Visualization

2. Geometrical Transformations

3. Setting Material Properties

You can use the SetEps and SetEpsMu class methods of RWGGeometry to override the MATERIAL designations in .scuffgeo files.

The prototypes are

    void RWGGeometry::SetEps(const char *Label, cdouble Eps, cdouble Mu);
    void RWGGeometry::SetEps(cdouble Eps);

    void RWGGeometry::SetEpsMu(const char *Label, cdouble Eps, cdouble Mu);
    void RWGGeometry::SetEpsMu(cdouble Eps, cdouble Mu);

The Label field selects the object you want to modify. (This is the label assigned using the OBJECT keyword in the .scuffgeo file.) The prototypes with no Label parameter are used to modify the exterior medium.

C++ example:

    RWGGeometry *G = new RWGGeometry("MyGeometry.scuffgeo");

    G->SetEps(3.4); // set the exterior medium to Eps = 3.4

    cdouble MyEpsValue(4.5,6.7);
    G->SetEps("OuterSphere", MyEpsValue); // set the object labeled OuterSphere to Eps = 4.5+6.7i


Note that, as soon as you make any call to SetEps or SetEpsMu, any user-specified or tabulated material designation you might previously have specified for the material in question is eliminated; the material properties remain constant until the next time you change them with SetEps or SetEpsMu, even if you subsequently reassemble the BEM matrix at a different frequency.

4. Specialized routines for computing field data

4A. Power, Force, and Torque Transferred to Objects by the incident field

Once we have solved a scattering problem, we can compute the power absorbed by a scattering object from the incident field and the force and torque exerted by the external field on the object.

The former quantity is obtained by integrating the Poynting vector over a surface bounding the object, while the latter quantities are obtained by integrating the Maxwell stress tensor over a surface bounding the object. Thus one way to get at these quantities would be to evaluate numerical cubatures over bounding surfaces, with the E and H fields at each point obtained by calling GetFields.

Instead, scuff-em suite takes advantage of certain specialized BEM tricks to compute these quantities much more directly and efficiently. The RWGGeometry class exports a routine named GetPFT that implements this calculation. (PFT stands for power, force, and torque.)

The prototypes are

   void RWGGeometry::GetPFT(HVector *KN, HVector *RHS, cdouble Omega, 
                            int SurfaceIndex, double PFT[8]);
   void RWGGeometry::GetPFT(HVector *KN, HVector *RHS, cdouble Omega, 
                            char *SurfaceLabel, double PFT[8]);

The input parameters are

  • KN: the solution vector of the BEM scattering problem.
  • RHS: the right-hand-side vector of the BEM scattering problem.
  • Omega: the angular frequency.
  • SurfaceIndex / SurfaceLabel: either the integer-valued index (first prototype) or the string-valued label (second prototype) of the surface (object) for which you want to compute the power, force, and torque.

    The PFT parameter must point to an array with room for at least 8 doubles. On return, this array is filled in as follows:

    Element Significance Unit
    PFT[0] Absorbed power Watts
    PFT[1] Total power Watts
    PFT[2] X force nanonewtons
    PFT[3] Y force nanonewtons
    PFT[4] Z force nanonewtons
    PFT[5] X torque nanonewtons • μm
    PFT[6] Y torque nanonewtons • μm
    PFT[7] Z torque nanonewtons • μm

    Note: Just to clarify, the element PFT[1] in the return vector is the total power, not the scattered power. Be careful, because this does not quite agree with what you might expect based on the way the PFT output files are written by scuff-scatter.

    Power scattered by objects

    The ``total power'' returned by GetPFT is the sum of the absorbed and scattered power. Hence the scattered power may be obtained as the difference between PFT[1] and PFT[0]. However, in cases where the total power is dominated by absorption, this may be an inaccurate way of computing the scattered power. An alternative method for calculating the scattered power, which does not suffer from the same inaccuracy, is provided by the GetScatteredPower routine.

       cdouble RWGGeometry::GetScatteredPower(HVector *KN, cdouble Omega, int SurfaceIndex);
       cdouble RWGGeometry::GetScatteredPower(HVector *KN, cdouble Omega, char *SurfaceLabel);

    The input parameters here have the same significance as their counterparts in GetPFT. The return value is the scattered power from the surface in question.

    Note that GetScatteredPower is more time-consuming than the method of GetPFT and hence should not be used unless you are in the aforementioned absorption-dominated regime.

    C++ example:

      /*- initialize the geometry ------------------------------------*/
      RWGGeometry *G=new RWGGeometry("MyGeometry.scuffgeo");
      /*- assemble and factorize the BEM matrix ----------------------*/
      double Omega=0.1;
      HMatrix *M   = G->AllocateBEMMatrix();
      G->AssembleBEMMatrix(Omega, M);
      /*- create a plane wave field source and assemble the RHS vector*/
      cdouble E0[3] = {1.0, 0.0, 0.0};
      double nHat[3] = {0.0, 0.0, 1.0};
      PlaneWave PW( E0, nHat );
      HVector *KN  = G->AllocateRHSVector();
      HVector *RHS = G->AllocateRHSVector();
      G->AssembleRHSVector(Omega, &PW, RHS);
      /*- solve the scattering problem     ---------------------------*/
      /*- compute power, force, torque on the object labeled "Cube"   */
      double PFT[8];
      G->GetPFT(KN, RHS, Omega, "Cube", PFT);
      printf("Absorbed power: %e watts \n", PFT[0]);
      printf("Total power:    %e watts \n", PFT[1]);
      printf("X force:        %e nanoNewtons \n", PFT[2]);
      printf("Z force:        %e nanoNewtons \n", PFT[4]);
      /*- compute scattered power ------------------------------------*/
      double PScat = PFT[1] - PFT[0];
      if ( fabs(PScat) < 0.01*fabs(PFT[1]) ) 
       { // PScat as computed via PFT methods is inaccurate;
         // recompute using more acccurate approach
         PScat = GetScatteredPower(KN, Omega, "Cube");
      printf("Scattered power: %e watts \n", PScat);

    4B. Dyadic Green's functions

    The dyadic Green's functions (DGFs) of a material configuration are certain susceptibilities that characterize the response of the configuration to illumination by point sources.

    DGFs may be computed by solving scattering problems in which we place a point source at some point x and compute the fields at some other point x. The total DGF is obtained by computing the total field at x. The scatttering DGF (or the "scattering part" of the DGF) is obtained by computing just the scattered field at x, neglecting the direct contribution of the original point source.

    More specifically, the scattering electric and magnetic DGFs are 3x3 matrices whose components have the following significance:

    The fact that BEM techniques naturally separate the contributions of induced sources from the contributions of incident-field sources makes the BEM particularly well-suited to the computation of scattering DGFs.

    A quantity of particular interest in many branches of physics is the scattering DGF evaluated for coincident source and evaluation points x=x. libscuff offers the built-in class method GetDyadicGFs() for computing this quantity at arbitrary coincident points in a scuff-em geometry. The prototype for this routine is this:

       void RWGGeometry::GetDyadicGFs(double X[3], cdouble Omega, HMatrix *M, HVector *KN, cdouble GE[3][3], cdouble GM[3][3]);

    The input parameters are

    • X: the Cartesian coordinates of the evaluation point
    • Omega: the angular frequency
    • M: The assembled and LU-factorized BEM matrix. It is the caller's responsibility to call AssembleBEMMatrix() and LUFactorize() before the call to GetDyadicGFs.
    • KN: A scratch workspace HVector used internally within the routine. KN should be the result of a previous call to AllocateRHSVector.

    On return, the GE and GM arrays are filled in with the components of the scattering electric and magnetic DGFs at the given evaluation point.

    5. Routines for evaluating numerical cubatures over panels

    In many cases it is convenient to evaluate surface integrals over the scattering surfaces in your geometry. For example, such integrals may be used to compute (1) the contribution of surface currents to the electromagnetic fields at arbitrary points in space, or (2) the electric and magnetic dipole moments induced by an incident field on the surfaces in your geometry.

    Of course, many such calculations are performed automatically for you behind the scenes by specialized libscuff routines; for example, item (1) is done for you by GetFields(), while item (2) is done for you by GetDipoleMoments(). However, in some cases you may need to evaluate your own types of surface integrals. For this purpose, libscuff offers several convenience routines, as discussed below.

    Since surfaces in scuff-em are specified as unions of flat triangular panels, surface integrals reduce to panel integrals. Broadly speaking, there are two cases of interest:

    1. You want to evaluate an integral over a single panel, i.e. you want to perform a panel cubature. This case arises, for example, in computing the multipole moments of the source distribution on an individual panel.
    2. You want to evaluate integrals describing the interaction between two panels, i.e. you want to perform a panel-panel cubature. For example, you may want to compute the average value of the electric field due to sources on one panel, with the average taken over the area of the other panel. This case arises in computing the elements of the matrices that enter into BEM solvers.

    More specifically, a panel cubature is a two-dimensional integral of the form

    where P is a panel, and where we have written the integrand to indicate that it may depend not only on the integration point x but also on the value b(x) of the RWG basis function at x. Similarly, a panel-panel cubature is a four-dimensional integral of the form

    where P1 and P2 are two panels, and where again the integrand may depend on the values of the RWG basis functions on the two panels.

    Closely related to the notion of a panel cubature is the notion of an integral over the support of an RWG basis function, which we term a basis-function cubature (or BF cubature for short). Since each RWG basis function is supported on two triangular panels, a basis-function cubature is just a sum of two panel cubatures. Similarly, an obvious extension of the panel-panel cubature is the BF-BF cubature, which is a four-dimensional integral over the supports of two RWG basis functions, computed in practice as a sum of four panel-panel integrals.

    libscuff offers the following routines for addressing each of the above situations.

        void GetPanelCubature(RWGSurface *S, int np, 
                              PCFunction Integrand, void *UserData, int IDim,
                              int Order, double RelTol, double AbsTol,
                              cdouble Omega, HVector *KN,
                              double *Result);
        void GetBFCubature(RWGSurface *S, int ne, 
                           PCFunction Integrand, void *UserData, int IDim,
                           int Order, double RelTol, double AbsTol,
                           cdouble Omega, HVector *KN,
                           double *Result);
        void GetPanelPanelCubature(RWGSurface *S1,  int np1, RWGSurface *S2, int np2,
                                   PPCFunction Integrand, void *UserData, int IDim,
                                   int Order, double RelTol, double AbsTol,
                                   cdouble Omega, HVector *KN,
                                   double *Result);
        void GetBFBFCubature(RWGSurface *S1,  int ne1, RWGSurface *S2, int ne2,
                             PPCFunction Integrand, void *UserData, int IDim,
                             int Order, double RelTol, double AbsTol,
                             cdouble Omega, HVector *KN,
                             double *Result);

    The various parameters here are explained in the table below.

    Parameter Meaning
    RWGSurface *S
    int np
    For the panel cubature case, these parameters specify the panel over which we integrate: S is the surface to which the panel belongs, and np is the index of the panel within that surface (0 <= np <= S->NumPanels-1).
    RWGSurface *S
    int ne
    For the BF cubature case, these parameters specify the RWG basis function over which we integrate: S is the surface to which the basis function belongs, and ne is the index of the mesh edge within that surface to which the basis function in question is associated (0 <= ne <= S->NumEdges-1).
    RWGSurface *S1
    int np1
    RWGSurface *S2
    int np2
    For the panel-panel cubature case, (S1,np1) and (S2,np2) identify the panels over which we integrate.

    For the BF-BF cubature case, (S1,ne1) and (S2,ne2) identify the mesh edges associated to the two basis functions over whose supports we integrate.

    PCFunction *Integrand User-specified integrand function for the panel cubature routine (see below).
    PPCFunction *Integrand User-specified integrand function for the panel-panel cubature routine (see below).
    void *UserData Optional user data passed to Integrand.
    int IDim The dimension of the integrand vector, measured in number of double values. Thus, if your Integrand routine computes one double value, set IDim=1; if your routine computes one cdouble value, set IDim=2; if your routine computes 6 different double-valued quantities, set IDim=6.
    int Order The order of the cubature scheme to use. Higher orders yield higher accuracy but are more computationally expensive. At present only two nonzero orders are supported: Order=4 and Order=20.

    Setting Order=0 will request an adaptive cubature that will recursively subdivide the integration domain until a desired error tolerance is achieved.

    double RelTol
    double AbsTol
    Relative and absolute error tolerances for the adaptive integration. These are only used if Order=0.
    cdouble Omega
    HVector *KN
    These optional parameters may be used to specify information about the surface-current distribution resulting from the solution to a scattering problem.

    If you have solved a scattering problem to obtain a vector of surface-current expansion coefficients KN, then you may pass that vector together with Omega into the panel cubature routines. In this case, the K,N,Eps,Mu fields in the PanelData structure passed to your integrand routine will be properly filled in.

    Alternatively, you may pass NULL for KN (and 0 or any arbitrary value for Omega), in which case your integrand routine just ignore the corresponding fields in the PanelData structure. In this case, there is no underlying surface-current distribution available, but the panel cubatures can still be used to compute geometric data on panels.

    double *Result Pointer to a user-allocated buffer with enough room to store IDim values of type double. On return, this vector is filled in with the values of the integral.

    Core Library

  • libscuff documentation: Ancillary RWGGeometry Operations, by Homer Reid
    Last Modified: 11/16/16