CGAL Mesh_3: How to adhere to a surface inside the domain using Polyhedral_complex_mesh_domain_3?

Question

The goal

My input is:

• several patches forming the boundary of the domain and
• one or more surface inside the domain.

What I would like to obtain is a tetrahedral mesh that:

1. fills the entire domain,
2. conforms to the interior surface and
3. its boundary vertices have an ID corresponding to their input patch.

One possibility would be to merge the input patches into one mesh and use the example Remeshing a polyhedral domain with surfaces from the documentation; this works well but fills only requirements 1 and 2.

Another possibility is, having learnt here how to keep the patch IDs, to modify the example to use Polyhedral_complex_mesh_domain_3 instead of Polyhedral_mesh_domain_with_features_3 so as to have the access to the constructor with subdomains. So far, my solution (see below) fills requirements 2 and 3 but not 1.

Code so far

As a simple example, I have split the file horizons-domain.off in two; the first file contains the first 10, the second the remaining 2 triangles.

sides.off

OFF
8 10 0

-1.1855500570497703 -0.076163891881438378 -0.8013403915768772
-1.1855500570497703 0.47597074519009164 -0.8013403915768772
0.79704321809070222 0.47597074519009164 -0.8013403915768772
0.79704321809070222 -0.076163891881438378 -0.8013403915768772
-1.1855500570497703 -0.076163891881438378 1.0953134363531141
-1.1855500570497703 0.47597074519009164 1.0953134363531141
0.79704321809070222 0.47597074519009164 1.0953134363531141
0.79704321809070222 -0.076163891881438378 1.0953134363531141
3  0 1 3
3  3 1 2
3  0 4 1
3  1 4 5
3  3 2 7
3  7 2 6
3  4 0 3
3  7 4 3
3  6 4 7
3  6 5 4

top.off

OFF
4 2 0

-1.1855500570497703 0.47597074519009164 -0.8013403915768772
0.79704321809070222 0.47597074519009164 -0.8013403915768772
-1.1855500570497703 0.47597074519009164 1.0953134363531141
0.79704321809070222 0.47597074519009164 1.0953134363531141
3  0 2 3
3  1 0 3

subdomains.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedral_complex_mesh_domain_3.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/make_mesh_3.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Polyhedral_complex_mesh_domain_3<K> Mesh_domain;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;

typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
    Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index> C3t3;

using namespace CGAL::parameters;

int main(int argc, char*argv[])
{
    // Read patches
    std::cout.precision(17);
    std::cerr.precision(17);

    std::ifstream input1(argv[1]);
    std::ifstream input2(argv[2]);
    std::ifstream input3(argv[3]);

    std::vector<Polyhedron> patches(3);
    input1 >> patches[0];
    input2 >> patches[1];
    input3 >> patches[2];

    // The first mesh is inside subdomain 0, the other two are on the boundary between subdomains 0 and 1.
    std::vector<std::pair<int, int>> incident_subdomains(3);
    incident_subdomains[0] = std::make_pair(0, 0);
    incident_subdomains[1] = std::make_pair(0, 1);
    incident_subdomains[2] = std::make_pair(0, 1);
    
    Mesh_domain domain(patches.begin(), patches.end(), incident_subdomains.begin(), incident_subdomains.end());

    // Mesh generation
    Mesh_criteria criteria(facet_distance=0.01, cell_radius_edge_ratio = 2);
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude());

    // Output
    dump_c3t3(c3t3, "out");
}

When I compile and run with

> ./subdomains horizons.off sides.off top.off

(horizons is the example file from the documentation), the tetrahedra try to conform to the interior surface, keep the IDs but do not fill the entire domain:

I also tried changing the line

incident_subdomains[0] = std::make_pair(0, 0);

to

incident_subdomains[0] = std::make_pair(1, 1);

This leads to a slightly better result but still far from requirement 3.

How can I fix it? I have tried various combinations of 1s and 0s in the std::pairs of the subdomains without success.