Triangulate faces while preserving custom data in CGAL

I have a Polyhedron that uses a custom Polyhedron_items_ 3 class to add custom data to the faces of the structure. For example:

    typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
    typedef Kernel::Point_3 Point_3;
    typedef Kernel::Vector_3 Vector_3;
    typedef Kernel::Plane_3 Plane_3;

    template <class Refs, class Traits>
    struct TFace : public CGAL::Polyhedron_items_3::Face_wrapper<Refs,Traits>::Face {
        std::array<double, 8> my_data;

        TFace() {}
        TFace(const Plane_3& p) : CGAL::Polyhedron_items_3::Face_wrapper<Refs,Traits>::Face(p) {}
    };

    struct CustomPolyhedronItems : public CGAL::Polyhedron_items_3 {
        template <class Refs, class Traits>
        struct Face_wrapper
        {
            typedef typename Traits::Plane_3 Plane;
            typedef TFace<Refs, Plane> Face;
        };
    };

    typedef CGAL::Polyhedron_3<Kernel, CustomPolyhedronItems> Polyhedron;
    typedef boost::graph_traits<Polyhedron>::vertex_descriptor vertex_descriptor;
    typedef boost::graph_traits<Polyhedron>::face_descriptor face_descriptor;

This works, and allows me to assign custom data (the my_data array above) to each face of a Polyhedron.

However, an operation like CGAL::Polygon_mesh_processing::triangulate_faces will generate new faces. Upon testing, it seems that the newly generated faces do not copy over the data from the face that was split, so my_data is not preserved.

To test, I created a simple cube (6 quads) with custom data on every face, ran it through triangulate_faces, and checked the custom data. Half of the triangles retained the original data, the other half were all zeroed.

Is it possible to tell triangulate_faces to copy my custom data when generating new faces? E.g. if a quad is converted into two triangles, then those triangles should have the same data as the original quad.

Also, if it is possible, would it also be possible to accomplish something similar with corefine_and_compute_union or something like a remesh operation?

EDIT:

In case that wasn't clear, here's a quick illustration of the problem using normals in place of my_data. Half of the triangles have correct normals, but the other half are wrong since they're zeroed:

Obviously, in the case of normals I can just recalculate them. But for my custom data, that's not possible. I need to make sure I preserve it when triangulating (and corefine union)

You have to use a Triangle Visitor. Unfortunately, I did that a couple of months ago and I don't remember everything. Note that my triangle mesh has not only face colors properties, but also face scalars.

struct TriangulateVisitor : 
  public PMP::Triangulate_faces::Default_visitor<EMesh3>
{
  void before_subface_creations(face_descriptor fsplit) {
    *ofaceindex = fsplit;
  }

  void after_subface_created(face_descriptor fnew) {
    (*fmap).insert(std::make_pair(fnew, *ofaceindex));
  }
  
  TriangulateVisitor()
    : fmap(new MapBetweenFaces()),
      ofaceindex(new face_descriptor())
  {}
  
  std::shared_ptr<MapBetweenFaces> fmap;
  std::shared_ptr<face_descriptor> ofaceindex;
};

where

typedef std::map<face_descriptor, face_descriptor> MapBetweenFaces;

and

void triangulateMesh(EMesh3& mesh) {
  MaybeFcolorMap fcolormap_ = 
    copy_prop<face_descriptor, std::string, EK>(mesh, "f:color");
  const bool hasFcolors = fcolormap_.second;
  MaybeFscalarMap fscalarmap_ = 
    copy_prop<face_descriptor, double, EK>(mesh, "f:scalar");
  const bool hasFscalars = fscalarmap_.second;
  removeProperties(mesh, {"v:normal"});
  TriangulateVisitor vis;
  const bool success = 
    PMP::triangulate_faces(mesh, CGAL::parameters::visitor(vis));
  if(!success) {
    Rcpp::stop("Triangulation has failed.");
  }
  if(hasFcolors) {
    MapBetweenFaces fmap = *(vis.fmap);
    Fcolors_map fcolors = 
      mesh.add_property_map<face_descriptor, std::string>(
        "f:color", ""
      ).first;
    for(EMesh3::Face_index fi: mesh.faces()) {
      fcolors[fi] = fcolormap_.first[fmap[fi]];
    }
  }
  if(hasFscalars) {
    MapBetweenFaces fmap = *(vis.fmap);
    Fscalars_map fscalars = 
      mesh.add_property_map<face_descriptor, double>(
        "f:scalar", nan("")
      ).first;
    for(EMesh3::Face_index fi: mesh.faces()) {
      fscalars[fi] = fscalarmap_.first[fmap[fi]];
    }
  }
}

Say me whether this is clear.