Using BoundaryIntegralOp class

The BoundaryIntegralOp class is designed to construct and evaluate boundary integral operators for solving problems in potential theory, such as the Laplace equation or the Stokes flow problem.

The following is a brief tutorial on using the BoundaryIntegralOp class. For more advanced usage and additional features, please refer to the class API in boundary_integral.hpp.

Building a Boundary Integral Operator

To build a boundary integral operator, follow these steps:

1. Initialize the Boundary Integral Operator

Initialize the boundary integral operator with the desired kernel and communication context. In this example, we use the Stokes double-layer kernel (Stokes3D_DxU).

const Stokes3D_DxU ker;
BoundaryIntegralOp<double,Stokes3D_DxU> BIOp(ker, false, comm);

2. Set Quadrature Accuracy

Set the quadrature accuracy for numerical integration.

BIOp.SetAccuracy(1e-10);

3. Add Geometry to the Operator

Add the geometry of the surface to the boundary integral operator. This includes specifying the discretization of the surface. In this example the boundary is given by an instance of SlenderElemList class (in CSBQ).

SlenderElemList<double> elem_lst;
elem_lst.Read<double>("data/loop.geom", comm);
BIOp.AddElemList(elem_lst);

4. Set Evaluation Points

The target points can be specified as follows. If not set or Xt is empty, then the default target points are the surface discretization nodes.

BIOp.SetTargetCoord(Xt);

Evaluating Potentials

Once the boundary integral operator is constructed, you can evaluate potentials at on- and off-surface target points.

1. Compute the Potential

Compute the potential using the boundary integral operator and a density function defined on the surface.

Vector<double> sigma(Ninput);
sigma = 1; // Set the density function sigma at each surface discretization node
Vector<double> U;
BIOp.ComputePotential(U, sigma); // Compute the potential U

2. Visualize the Results

Visualize the geometry and the computed potential.

elem_lst.WriteVTK("vis/U", U, comm); // Write visualization data to VTK file

Example Code

Below is an example code demonstrating the usage of the BoundaryIntegralOp class:

#include <sctl.hpp>
using namespace sctl;

int main(int argc, char** argv) {
  Comm::MPI_Init(&argc, &argv);

  {
    const Comm comm = Comm::World();

    const Stokes3D_DxU ker;
    BoundaryIntegralOp<double,Stokes3D_DxU> BIOp(ker, false, comm);
    BIOp.SetAccuracy(1e-10);

    SlenderElemList<double> elem_lst;
    elem_lst.Read<double>("data/loop.geom", comm); // load geometry
    BIOp.AddElemList(elem_lst); // add element list to boundary integral operator

    //BIOp.SetTargetCoord(Xt); // set target points (default is discretization nodes)

    const Long Ninput = BIOp.Dim(0); // (local) input dimension of the operator

    Vector<double> sigma(Ninput);
    sigma = 1;

    Vector<double> U;
    BIOp.ComputePotential(U, sigma); // compute potential

    elem_lst.WriteVTK("vis/U", U, comm); // write visualization
  }

  Comm::MPI_Finalize();
  return 0;
}