BRepCheck_SubshapeNotInShape error

[FlachyJoe]

Background:
I'm writing an algorithm to deform a shape.
for each face

1. emptyCopy the initial face
2. for each wire of face
3. discretize edges
4. move the points
5. create a BSpline curve on them
6. BRepBuilderAPI_MakeEdge
7. BRepBuilderAPI_MakeWire the new wire
8. add the wires on the new face with TopoDS_Builder::Add
9. ShapeFix_Wire

then
TopoDS_Builder::MakeShell an empty shell
TopoDS_Builder::Add the faces in the shell
BRepBuilderAPI_Sewing the shell
BRepBuilderAPI_MakeSolid
ShapeFix_Solid

Issue:
BRepCheck on the result returns BRepCheck_SubshapeNotInShape

Why and how can I fix this?

[dpasukhi]

Please provide a code sample and result shape (.brep file, shared as zip or renamed to .txt).

[FlachyJoe]

Input and output of the code:
input-cube.brep.txt
output-cube.brep.txt

/*
 * (c) 2026 Florian Foinant-Willig <ffw_at_2f2v_dot_fr>
 */

#include <iostream>
#include <sys/stat.h>

// OpenCascade includes
#include <Approx_Curve3d.hxx>

#include <BRep_Builder.hxx>

#include <BRepAdaptor_Curve.hxx>

#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepBuilderAPI_MakeSolid.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_Sewing.hxx>

#include <BRepCheck_Analyzer.hxx>

#include <BRepTools.hxx>
#include <BRepTools_ReShape.hxx>

#include <GCPnts_UniformAbscissa.hxx>

#include <Geom_BSplineCurve.hxx>
#include <Geom_TrimmedCurve.hxx>

#include <GeomAPI_Interpolate.hxx>
#include <GeomAPI_PointsToBSpline.hxx>

#include <gp_Pnt.hxx>

#include <ShapeFix_Shape.hxx>
#include <ShapeFix_Wire.hxx>

#include <TopoDS.hxx>
#include <TopoDS_Builder.hxx>

#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>


using namespace std;

// check file existence
inline bool fileExists (const std::string& name) {
    struct stat buffer;
    return (stat (name.c_str(), &buffer) == 0);
}

// Read BRep from file
TopoDS_Shape importBREP(string filename){
    TopoDS_Shape s;
    BRep_Builder b;
    std::ifstream is;
    is.open(filename);
    BRepTools::Read(s, is, b);
    is.close();
    return s;
}

// Write BRep to file
void exportBREP(string filename, TopoDS_Shape s){
    BRep_Builder b;
    std::ofstream os;
    os.open(filename);
    BRepTools::Write(s, os);
    os.close();
}

// apply a deformation to a point
gp_Pnt tensorTransform(gp_Pnt from){
    auto pi_40 = atan(1.)/10;
    auto x = from.X() ;
    auto y = from.Y() * cos(pi_40 * from.X()) - from.Z() * sin(pi_40 * from.X());
    auto z = from.Z() * sin(pi_40 * from.X()) + from.Z() * cos(pi_40 * from.X());

    return gp_Pnt(x, y, z);
}


// from FreeCAD GeomBSplineCurve::approximate
Handle(Geom_BSplineCurve) approximate(TColgp_Array1OfPnt pnts, double tol3d, int minDegree, int maxDegree, GeomAbs_Shape continuity)
{
    GeomAPI_PointsToBSpline fit(pnts, Approx_ParametrizationType::Approx_ChordLength, minDegree, maxDegree, continuity, tol3d);
    return fit.Curve();
}

// initialize a simple curve
Handle(Geom_BSplineCurve) newBSplineCurve() {
    TColgp_Array1OfPnt poles(1, 2);
    poles(1) = gp_Pnt(0.0, 0.0, 0.0);
    poles(2) = gp_Pnt(1.0, 0.0, 0.0);
    TColStd_Array1OfReal knots(1, 2);
    knots(1) = 0.0;
    knots(2) = 1.0;
    TColStd_Array1OfInteger mults(1, 2);
    mults(1) = 2;
    mults(2) = 2;
    return new Geom_BSplineCurve(poles, knots, mults, 1);
}


// tweak of FreeCAD TopoShapeFacePy::validate
void validateFace(TopoDS_Face& face)
{
    BRepCheck_Analyzer aChecker(face);
    if (!aChecker.IsValid()) {
        TopoDS_Wire outerwire = BRepTools::OuterWire(face);

        TopExp_Explorer xp(face, TopAbs_WIRE);
        ShapeFix_Wire fix;
        fix.SetFace(face);
        fix.Load(outerwire);
        fix.ModifyTopologyMode() = Standard_False;
        fix.ModifyGeometryMode() = Standard_True;
        fix.Perform();
        BRepBuilderAPI_MakeFace mkFace(fix.Wire());
        while (xp.More()) {
            if (outerwire.IsNotEqual(xp.Current())) {
                fix.Load(TopoDS::Wire(xp.Current()));
                fix.ModifyTopologyMode() = Standard_False;
                fix.ModifyGeometryMode() = Standard_True;
                fix.Perform();
                mkFace.Add(fix.Wire());
            }
            xp.Next();
        }

        aChecker.Init(mkFace.Face());
        if (!aChecker.IsValid()) {
            ShapeFix_Shape fix(mkFace.Face());
            fix.SetPrecision(Precision::Confusion());
            fix.SetMaxTolerance(Precision::Confusion());
            fix.Perform();
            fix.FixWireTool()->Perform();
            fix.FixFaceTool()->Perform();
            auto result = fix.Shape();
            if (result.ShapeType() == TopAbs_FACE) {
                face = TopoDS::Face(fix.Shape());
                return;
            }
        }
        face = TopoDS::Face(mkFace.Face());
    }
}

TopoDS_Wire deformWire(TopoDS_Wire wire) {
    BRepBuilderAPI_MakeWire wireBuilder;

    vector<TopoDS_Edge> newEdges;

    // for each edge
    for(TopExp_Explorer exp(wire, TopAbs_EDGE); exp.More(); exp.Next())
    {
        auto edge = TopoDS::Edge(exp.Current());

        // discretize
        BRepAdaptor_Curve adapt(edge);
        int numPoints = 10;
        double first = adapt.FirstParameter();
        double last = adapt.LastParameter();
        GCPnts_UniformAbscissa discretizer;
        discretizer.Initialize(adapt, numPoints, first, last);
        int nbPoints = discretizer.NbPoints();
        TColgp_Array1OfPnt points(1, nbPoints);
        TColgp_Array1OfPnt deformPoints(1, nbPoints);

        //deform
        for (int i = 1; i <= nbPoints; i++) {
            auto p = adapt.Value(discretizer.Parameter(i));
            points(i) = p;
            deformPoints(i) = tensorTransform(p);
        }

        // approximate the points
        auto resultCurve = approximate(deformPoints, 1e-6, 3, 8, GeomAbs_Shape::GeomAbs_C2);
        if (adapt.IsPeriodic() || adapt.IsClosed()) {
            resultCurve->SetPeriodic();
        }

        // recreate an edge
        BRepBuilderAPI_MakeEdge mkBuilder(resultCurve, resultCurve->FirstParameter(), resultCurve->LastParameter());
        auto ne = TopoDS::Edge(mkBuilder.Shape());
        ne.Orientation(edge.Orientation());
        newEdges.push_back(ne);
    } // edge

    for (auto& ne : newEdges) {
        wireBuilder.Add(ne);
    }
    wireBuilder.Build();
    if (wireBuilder.IsDone()) {
        return wireBuilder.Wire();
    }

    auto newWire = TopoDS::Wire(wire.EmptyCopied());
    TopoDS_Builder builder;
    for (auto& ne : newEdges) {
        builder.Add(newWire, ne);
    }
    return newWire;
}

int main(int argc, char *argv[])
{
    if (argc != 3) {
        cout << "Usage :\n"
             << argv[0] << " input.brep output.brep\n" ;
        return -1;
    }

    auto inputFile = argv[1];
    if (! fileExists(inputFile)) {
        cout << argv[1] << " not found.\n";
        return -1;
    }

    auto shape = importBREP(inputFile);
    if (shape.IsNull()){
        cout << "Null input shape.\n";
        return -1;
    }

    TopoDS_Builder topoBuilder;

    TopoDS_Shell resultShell;
    topoBuilder.MakeShell(resultShell);

    //for each face
    for(TopExp_Explorer expF(shape, TopAbs_FACE); expF.More(); expF.Next())
    {
        auto face = TopoDS::Face(expF.Current());
        auto resultFace = TopoDS::Face(face.EmptyCopied());

        auto outerwire = BRepTools::OuterWire(face);
        topoBuilder.Add(resultFace, deformWire(outerwire));

        for(TopExp_Explorer exp(face, TopAbs_WIRE); exp.More(); exp.Next())
        {
            auto wire = TopoDS::Wire(exp.Current());

            if (wire.IsNotEqual(outerwire)) {
                auto w = deformWire(wire);
                topoBuilder.Add(resultFace, w);
            }
        }
        validateFace(resultFace);
        topoBuilder.Add(resultShell, resultFace);
    }


    BRepBuilderAPI_MakeSolid solidifer;
    solidifer.Add(resultShell);

    ShapeFix_Solid fixer(solidifer.Solid());
    auto resultSolid = fixer.Solid();

    exportBREP(argv[2], resultSolid);

    // BRepBuilderAPI_Sewing sewer;
    // sewer.Add(resultShell);
    // auto sewed = sewer.SewedShape();
    // if (sewed.IsNull()) {
    //     if (resultShell.IsNull()) {
    //         cout << "Result shape is null !\n";
    //         return -1;
    //     }
    //     cout << "Failed to sew\n";
    //     exportBREP(argv[2], resultShell);
    // }
    // else {
    //     exportBREP(argv[2], sewed);
    // }

    cout << "Done!\n";
    return 0;
}

[FlachyJoe]

This problem was due to a lack of sewing of the shell before the solidification.