Finding all connected pieces from an implicit function

[yangxucmt]

Dear all,

I am quite new to CGAL but I think this is a common question so was hoping it could be solved very easily.

I am trying to generate a triangulation of a 2d surface determined by the implicit function f(x,y,z)=E. The issue with my system is that depending on the value of E, the surface is either one connected piece or two connected piece. And the code attached can only find one connected piece.

Although I have some analytic understanding of the two components and when it appears, I would like my code to automatically find all the connected pieces and output all the surface mesh triangles that can be further utilized. My code was below and I know it is not the most efficient, so any advice is welcome.

Best,
Xu

Code:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Labeled_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>

#include <CGAL/IO/File_medit.h>

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <set>

namespace params = CGAL::parameters;

//======================================================
// Kernel
//======================================================

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = K::FT;
using Point = K::Point_3;

//======================================================
// Domain
//======================================================

using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;

#ifdef CGAL_CONCURRENT_MESH_3
using Concurrency_tag = CGAL::Parallel_tag;
#else
using Concurrency_tag = CGAL::Sequential_tag;
#endif

using Tr =
    CGAL::Mesh_triangulation_3<
        Mesh_domain,
        CGAL::Default,
        Concurrency_tag>::type;

using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;

using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;

using Facet_criteria = Mesh_criteria::Facet_criteria;
using Cell_criteria = Mesh_criteria::Cell_criteria;

//======================================================
// Physical parameters
//======================================================

const double m = 0.2;
const double lambda = 4.0;

const double vx = 10.0 * std::sqrt(lambda) / m;
const double vz = 10.0 * std::sqrt(lambda) / m;
const double ux = 0.8 * vx;

const double muvalue = 6.0;

//======================================================
// Implicit function
//======================================================

FT f(const Point &p)
{
    const FT kx = p.x();
    const FT ky = p.y();
    const FT kz = p.z();

    const FT ky2 = ky * ky;

    FT energy =
        (2.0 * m * ux * kx +
         std::sqrt(
             ky2 * ky2 - 2.0 * lambda * ky2 + lambda * lambda + 4.0 * m * m * vx * vx * kx * kx + 4.0 * m * m * vz * vz * kz * kz)) /
        (2.0 * m);

    return energy - muvalue;
}

//======================================================

int main()
{
    //--------------------------------------------------
    // Domain
    //--------------------------------------------------

    Mesh_domain domain =
        Mesh_domain::create_implicit_mesh_domain(
            f,
            K::Sphere_3(CGAL::ORIGIN, 100.0),
            params::relative_error_bound(1e-8));

    //--------------------------------------------------
    // Criteria
    //--------------------------------------------------

    Facet_criteria facet_criteria(
        30.0, // angle
        0.05, // size
        0.01 // approximation
    );

    Cell_criteria cell_criteria(
        2.0,
        0.1);

    Mesh_criteria criteria(
        facet_criteria,
        cell_criteria);

    //--------------------------------------------------
    // Mesh generation
    //--------------------------------------------------

    C3t3 c3t3 =
        CGAL::make_mesh_3<C3t3>(
            domain,
            criteria,
            params::surface_only()
                .no_exude()
                .no_perturb());

    //--------------------------------------------------
    // Save mesh
    //--------------------------------------------------

    std::ofstream medit("surface.mesh");
    CGAL::IO::write_MEDIT(medit, c3t3);

    std::ofstream off("surface.off");
    c3t3.output_boundary_to_off(off);

    //--------------------------------------------------
    // Statistics
    //--------------------------------------------------

    const Tr &tr = c3t3.triangulation();

    std::cout << std::setprecision(16);

    std::cout
        << "Vertices : "
        << tr.number_of_vertices()
        << std::endl;

    std::cout
        << "Surface facets : "
        << c3t3.number_of_facets_in_complex()
        << std::endl;

    std::cout
        << "Cells : "
        << c3t3.number_of_cells_in_complex()
        << std::endl;

    //--------------------------------------------------
    // Area
    //--------------------------------------------------
    //--------------------------------------------------
    // Area
    //--------------------------------------------------

    double area = 0.0;

    for (C3t3::Facets_in_complex_iterator
             fit = c3t3.facets_in_complex_begin();
         fit != c3t3.facets_in_complex_end();
         ++fit)
    {
        auto cell = fit->first;
        int id = fit->second;

        K::Point_3 p[3];

        int k = 0;
        for (int i = 0; i < 4; ++i)
        {
            if (i != id)
            {
                // Mesh_3 stores Weighted_point_3
                p[k++] = cell->vertex(i)->point().point();
            }
        }

        auto cross = CGAL::cross_product(p[1] - p[0],
                                         p[2] - p[0]);

        area += 0.5 * std::sqrt(cross.squared_length());
    }

    std::cout << std::setprecision(16);
    std::cout << "Area = " << area << std::endl;

    return 0;
}