Demo 2: Boundary Integral Operator

The code (tutorial/demo2-bio.cpp) demonstrates the basic usage of the BoundaryIntegralOp class. The BoundaryIntegralOp class (in SCTL) is designed to construct and evaluate boundary integral operators, such as the Stokes single-layer and double-layer potential operators. 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 with geometry data loaded from a file.

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

Compiling and Running the Code

Without MPI: Navigate to the project root directory and run:
make bin/demo2-bio && ./bin/demo2-bio
With MPI: In the project root directory, edit Makefile to set CXX=mpicxx and uncomment CXXFLAGS += -DSCTL_HAVE_MPI. Then, run:
make bin/demo2-bio && mpirun -n 2 --map-by slot:pe=1 ./bin/demo2-bio

This will evaluate the Stokes double-layer potential on a loop geometry and write the VTK visualization to the file vis/U.pvtu which can be opened in ParaView.

Complete Example Code

/**
 * This demo code shows how to use the class sctl::BoundaryIntegralOp. The
 * class is declared in SCTL/include/sctl/boundary_integral.hpp. In particular,
 * we show how to:
 *
 * 1) build a boundary integral operator.
 *
 * 2) evaluate potentials at on- and off-surface target points.
 *
 * To compile and run the code, start in the project root directory and run:
 * make bin/demo2-bio && ./bin/demo2-bio
 */

#include <csbq.hpp>
using namespace sctl;

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

  {
    const Comm comm = Comm::World(); // comm must have local scope since it needs to be destroyed before Comm::MPI_Finalize()

    Stokes3D_DxU ker; // Stokes double-layer kernel (@see SCTL/include/sctl/kernel_function.hpp)
    BoundaryIntegralOp<double,Stokes3D_DxU> BIOp(ker, false, comm); // construct the boundary integral operator
    BIOp.SetAccuracy(1e-10); // set quadrature accuracy

    SlenderElemList<double> elem_lst;
    elem_lst.Read<double>("data/ring.geom", comm); // read geometry data from file
    BIOp.AddElemList(elem_lst); // add geometry to the boundary integral operator

    // 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);

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

    // Set the density function sigma at each surface discretization node.
    // The nodes can be queried using elem_lst.GetNodeCoord() (@see ex1_geometry.cpp).
    Vector<double> sigma(Ninput);
    sigma = 1;

    Vector<double> U;
    BIOp.ComputePotential(U, sigma); // compute the potential U
    SCTL_ASSERT(U.Dim() == Noutput); // U will be of size Noutput

    // visualizing the geometry and the potential
    elem_lst.WriteVTK("vis/U", U, comm);
  }

  Comm::MPI_Finalize();
  return 0;
}