How to process output: trim_with_solid

trim_with_solid is called:

// Trim mesh A (m_VA and m_FA) with solid B (m_VB and m_FB)
igl::copyleft::cgal::trim_with_solid(
                                     m_VA, m_FA           // input: mesh A
                                     , m_VB, m_FB         // input: solid B
                                     , m_V, m_F, m_D, m_J // output
                                     );

The output is processed:

    // Loop over output vertex data
    for (int i = 0; i < m_V.rows(); ++i) {
        // Access all the vertex data
    }

   // Loop over output triangle (face) data
    for (int i = 0; i < m_F.rows(); ++i) {
        // Condition to decide if facet is inside or outside the solid B
        if (m_D.coeff(i)) {
            // The output face is outside the solid B
        } else {
            // The output face is inside the solid B
        }
    }

Fine result

Sometimes the trim result is fine. Like this one which is trimming a teapot model with a solid cylinder:

Bad result

But sometimes the trim result is NOT as expected:

As shown by the wire-frame view above, the intersection of teapot and cylinder is identified correctly. But for some reason, this loop is NOT able to detect triangles which are outside the cylinder:

   // Loop over output triangle (face) data
    for (int i = 0; i < m_F.rows(); ++i) {
        // Condition to decide if facet is inside or outside the solid B
        if (m_D.coeff(i)) {
            // Face is outside the solid B
        } else {
            // Face is inside the solid B
        }
    }

Does anybody know what I'm possibly missing?

Previously, I was using Qt3DExtras::QCylinderMesh as solid B (VB and FB). But it was NOT watertight. I replaced it with the following code which creates watertight cylinders.

Header:

    // Create a cylinder mesh which is exactly solid
    // without any open boundary
    // Down facing cylinder in Z direction with arguments:
    // r: radius
    // h: Height
    // ssteps: how many edges will create the base circle
    // sp: starting point
    static void cylinder(
            std::vector<QVector3D>& points
            , std::vector<std::array<int, 3>>& indices
            , const float r
            , const float h
            , const size_t ssteps = 45
            , const QVector3D &sp = QVector3D()
            );

Source:

void EditorUtils::cylinder(
        std::vector<QVector3D> &points
        , std::vector<std::array<int, 3>> &indices
        , const float r
        , const float h
        , const size_t ssteps
        , const QVector3D &sp
        )
{
    assert(ssteps > 0);

    auto PI = M_PI;

    auto steps = int(ssteps);
    points.reserve(2*ssteps);
    double a = 2*PI/steps;

    QVector3D jp = sp;
    QVector3D endp = {sp.x(), sp.y(), sp.z() + h};

    // Upper circle points
    for(int i = 0; i < steps; ++i) {
        double phi = i*a;
        double ex = endp.x() + r*std::cos(phi);
        double ey = endp.y() + r*std::sin(phi);
        points.emplace_back(ex, ey, endp.z());
    }

    // Lower circle points
    for(int i = 0; i < steps; ++i) {
        double phi = i*a;
        double x = jp.x() + r*std::cos(phi);
        double y = jp.y() + r*std::sin(phi);
        points.emplace_back(x, y, jp.z());
    }

    // Now create long triangles connecting upper and lower circles
    indices.reserve(2*ssteps);
    auto offs = steps;
    for(int i = 0; i < steps - 1; ++i) {
        indices.emplace_back(std::array<int, 3>{i, i + offs, offs + i + 1});
        indices.emplace_back(std::array<int, 3>{i, offs + i + 1, i + 1});
    }

    // Last triangle connecting the first and last vertices
    auto last = steps - 1;
    indices.emplace_back(std::array<int, 3>{0, last, offs});
    indices.emplace_back(std::array<int, 3>{last, offs + last, offs});

    // Try generating a watertight body.
    // So we create a triangle fan for the upper and lower
    // ending of the cylinder to close the geometry.
    points.emplace_back(jp); int ci = int(points.size() - 1);
    for(int i = 0; i < steps - 1; ++i)
        indices.emplace_back(std::array<int, 3>{i + offs + 1, i + offs, ci});

    indices.emplace_back(std::array<int, 3>{offs, steps + offs - 1, ci});

    points.emplace_back(endp); ci = int(points.size() - 1);
    for(int i = 0; i < steps - 1; ++i)
        indices.emplace_back(std::array<int, 3>{ci, i, i + 1});

    indices.emplace_back(std::array<int, 3>{steps - 1, 0, ci});

    return;
}

Now the holes are created properly: