How to create a periodic Isosurface in a periodic 3D Delaunay mesh?

[gy114]

Hi. I'm trying to find an isosurface defined by some known function F(x,y,z)=F0 in a 3D periodic Delaunay mesh.
Besides the built-in algorithms to generate a 3D periodic mesh, is there a straightforward way to implement this?
I'm thinking about marching through all the "unique cells". When any intersection is detected, I can interpolate one or two triangle pieces and collect them. But I realized several difficulties:

(1) It may not be straightforward to go through the edges of unique cells without any repeats. Can I use the iterator of unique edges? I realized that some of the edges in unique cells are not unique.

(2) I'm trying to calculate the Gaussian curvature of each vertex on the isosurface. That means I need access to the 1-ring. It seems that there's not a built-in data structure to describe a periodic surface mesh. Am I missing something?

[MaelRL]

You can generate a periodic surface mesh using the 3D mesher and not give any refinement criterion for cells.

(1) For periodic triangulations, you have access to unique simplices, see "Unique_..." in (for example) https://doc.cgal.org/latest/Periodic_3_triangulation_3/classCGAL_1_1Periodic__3__triangulation__3.html.

(2) You would use the 3D periodic triangulation as you would use a normal triangulation: you walk through simplices using the same cell->vertex(index) or cell->neighbor(index) etc. The surface will be defined by the subdomain indices of the cells of the mesh. You can have a look at https://doc.cgal.org/latest/SMDS_3/index.html#title0.

If you are just interested in extracting your periodic surface, have a look at https://github.com/CGAL/cgal/blob/master/SMDS_3/include/CGAL/facets_in_complex_3_to_triangle_mesh.h, and specifically facets_in_complex_3_to_triangle_mesh().

[gy114]

Hi, MaelRL.

I wrote the following code using the idea of facets_in_complex_3_to_triangle_mesh().
Then I followed your comment in #7932 trying to isotropically remesh the surface after scaling it. Now I'm testing the function without scaling (scaled by an identity arma::eye(3,3). The output is written in both vtp and off formats.
There's seems an issue. The figure shows the surface mesh after remeshing. Please help.

#include <vtkSmartPointer.h>
#include <vtkUnstructuredGrid.h>
#include <vtkDoubleArray.h>
#include <vtkPoints.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkPolygon.h>
#include <vtkXMLPolyDataWriter.h>


#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_triangulation_ds_vertex_base_3.h>
#include <CGAL/Periodic_3_triangulation_ds_cell_base_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Triangulation_cell_base_with_info_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_data_structure_3.h>
#include <CGAL/make_periodic_3_mesh_3.h>
#include <CGAL/Periodic_3_mesh_triangulation_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/facets_in_complex_3_to_triangle_mesh.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/surface_Delaunay_remeshing.h>



#include <CGAL/number_type_config.h> // CGAL_PI

#include <cmath>
#include <random>
#include <iostream>
#include <unordered_set>
#include <armadillo>




typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K> P3_Traits;
typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<> P3_VbDS;
typedef CGAL::Triangulation_vertex_base_3<P3_Traits, P3_VbDS >       P3_Vb;
typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<> P3_CbDS;
typedef CGAL::Triangulation_cell_base_3<P3_Traits, P3_CbDS >       P3_Cb;


typedef CGAL::Triangulation_vertex_base_with_info_3<double, P3_Traits, P3_Vb > VbInfo; // Test to attach data to vertices.
typedef CGAL::Triangulation_cell_base_with_info_3<double, P3_Traits, P3_Cb > CbInfo; // Test to attach data to cells.
typedef CGAL::Triangulation_data_structure_3<VbInfo, CbInfo >                 TDS;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<P3_Traits, TDS >               P3_Delaunay;
typedef P3_Delaunay::Point P3_Point;
typedef P3_Delaunay::Iso_cuboid Iso_cuboid;
typedef CGAL::Surface_mesh<P3_Point> Surface_mesh;


typedef CGAL::Labeled_mesh_domain_3<K>                      Periodic_mesh_domain;
// Triangulation
typedef CGAL::Periodic_3_mesh_triangulation_3<Periodic_mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr>                       C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr>                           Periodic_mesh_criteria;
// To avoid verbose function and named parameters call
namespace params = CGAL::parameters;



void write_to_vtp_SurfaceMesh(const std::string& filename, const Surface_mesh& mesh) {
    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
    vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();

    // Transfer vertices to VTK points
    for (auto v : mesh.vertices()) {
        const P3_Point& p = mesh.point(v);
        points->InsertNextPoint(p.x(), p.y(), p.z());
    }

    // Transfer faces to VTK polygons
    for (auto f : mesh.faces()) {
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        polygon->GetPointIds()->SetNumberOfIds(mesh.degree(f));
        // Add vertex indices to the polygon
        size_t i = 0;
        for (auto v : vertices_around_face(mesh.halfedge(f), mesh)) {
            polygon->GetPointIds()->SetId(i++, v.idx());
        }
        polygons->InsertNextCell(polygon);
    }

    polyData->SetPoints(points);
    polyData->SetPolys(polygons);

    // Write to VTK file
    vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
    writer->SetFileName(filename.c_str());
    writer->SetInputData(polyData);
    //writer->SetDataModeToBinary();
    writer->SetDataModeToAscii();
    writer->Write();
}

double periodic_function(const P3_Point& p) {
    return -std::cos(2 * CGAL_PI * p[0]) - std::cos(2 * CGAL_PI * p[1]) - std::cos(2 * CGAL_PI * p[2]);
}


struct halfedge2edge
{
    halfedge2edge(const Surface_mesh& m, std::vector<boost::graph_traits<Surface_mesh>::edge_descriptor>& edges)
            : m_mesh(m), m_edges(edges)
    {}
    void operator()(const boost::graph_traits<Surface_mesh>::halfedge_descriptor& h) const
    {
        m_edges.push_back(edge(h, m_mesh));
    }
    const Surface_mesh& m_mesh;
    std::vector<boost::graph_traits<Surface_mesh>::edge_descriptor>& m_edges;
};


double average_edge_length(const Surface_mesh& mesh) {
    double total_length = 0.0;
    unsigned int edge_count = 0;

    for (const auto& e : mesh.edges()) {
        auto h = mesh.halfedge(e);
        P3_Point p1 = mesh.point(mesh.source(h));
        P3_Point p2 = mesh.point(mesh.target(h));
        total_length += std::sqrt(CGAL::squared_distance(p1, p2));
        edge_count++;
    }

    return edge_count != 0 ? total_length / edge_count : 0.0;
}


int testStandardPeriodicMeshGenerationAndWriteToVTP_V2() {
    // Define the cube that bounds the periodic domain
    Iso_cuboid domain_bounds(P3_Point(-0.5, -0.5, -0.5), P3_Point(0.5, 0.5, 0.5));

    // Create domain using the periodic function directly
    Periodic_mesh_domain domain = Periodic_mesh_domain::create_implicit_mesh_domain(periodic_function, domain_bounds);
    Periodic_mesh_criteria criteria(params::facet_angle(30)
                                            .facet_size(0.05 )
                                            .facet_distance(0.025 )
                                            .cell_radius_edge_ratio(2.)
                                            .cell_size(0.05));
    C3t3 c3t3_bis = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria,  params::no_odt().
            no_lloyd().
            no_perturb().
            no_exude());
    // Output
    std::ofstream outfile1("outfile1.off");
    c3t3_bis.output_facets_in_complex_to_off(outfile1);
    std::cout << "Mesh has " << c3t3_bis.number_of_cells() << " cells." << std::endl;
    outfile1.close();
    Surface_mesh tmesh;
    CGAL::facets_in_complex_3_to_triangle_mesh(c3t3_bis, tmesh);

    //We now scale the old points to the non-orthogonal unit cell.
    /*arma::mat B = {{1.5064,  -0.8697,   0},
                   {1.5064,   0.8697,   0},
                   {0,        0,   0.9398}};*/

    arma::mat B = arma::eye(3,3);

    for (auto v : tmesh.vertices())  {
        P3_Point& p = tmesh.point(v);
        arma::vec armaP = {p.x(), p.y(), p.z()};
        arma::vec scaledArmaP = (B.t())*armaP;
        p = P3_Point(scaledArmaP(0), scaledArmaP(1), scaledArmaP(2));
    }
    std::cout << "Writting scaled surface to file... ";
    write_to_vtp_SurfaceMesh("outputSurfaceAfterScaling.vtp", tmesh);
    std::cout << "done." << std::endl;

    // Create a property map for edge constraints
    typedef boost::graph_traits<Surface_mesh>::edge_descriptor EdgeDescriptor;
    std::unordered_map<EdgeDescriptor, bool> constrained_edges;

    // Mark border edges as constrained
    int edgeCount = 0;
    for(EdgeDescriptor e : edges(tmesh)) {
        if(tmesh.is_border(e)) {
            constrained_edges[e] = true;
            edgeCount++;
        }
    }
    std::cout << edgeCount << " number of boarder edges are found.\n";



    double target_edge_length = average_edge_length(tmesh);
    std::cout << target_edge_length << " is the target edge length\n";
    // Apply isotropic remeshing with constraints on border edges
    CGAL::Polygon_mesh_processing::isotropic_remeshing(
            faces(tmesh),
            target_edge_length,
            tmesh,
            CGAL::parameters::edge_is_constrained_map(
                    boost::make_assoc_property_map(constrained_edges)
            ).protect_constraints(true)
            // other parameters...
    );
    std::cout << "Writting remeshed surface to file... ";
    write_to_vtp_SurfaceMesh("outputSurfaceAfterRemeshing.vtp", tmesh);
    std::cout << "done." << std::endl;

    std::ofstream off_file("meshAfterRemeshing.off");
    CGAL::write_off(off_file, tmesh);
    off_file.close();

    return 0;
}

int main() {
    testStandardPeriodicMeshGenerationAndWriteToVTP_V2();
    return 0;
}