how to save a periodical mesh.

[luogyong]

I used the example project geogram_demo_Delaunay3d (with appropriate modifications) to generate a periodic mesh and then output it in the mesh file format. However, I found two issues with the output file:

1) The coordinates of the output nodes are very strange, extremely large, far from the region, and there are duplicates, such as:

Vertices
409
...
-6.27743856e+66 -6.27743856e+66 -6.27743856e+66 0 
-6.27743856e+66 -6.27743856e+66 -6.27743856e+66 0 
-6.27743856e+66 -9.74583694e+297 -1.45681599e+144 0 
2.98671411e-307 1.28494769e-311 1.28504795e-311 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
-1.45681599e+144 -1.45681599e+144 -1.45681599e+144 0 
3.45813628e-230 1.28504859e-311 1.28505025e-311 0 
...

2) Many of the tetrahedral unit node numbers in the output are greater than the number of nodes (nb_vertices) output above, such as:

Tetrahedra
1002
10 91 12 87 0 
49 31 53 50 0 
71 28 78 44 0 
70 9 80 4 0 
1599 (>409) 1533 (>409) 7 1622 (>409) 0 
...

My question is, how to correctly input a periodic mesh? The code I used in the save function is as follows:

bool DemoDelaunay3dApplication::save(const std::string& filename) {
    MeshIOFlags flags;

    std::vector<double> pts(delaunay_->nb_vertices() * 3);
    std::vector<index_t> tri2v;
    for (index_t v = 0; v < delaunay_->nb_vertices(); ++v) {
        pts[3 * v] = delaunay_->vertex_ptr(v)[0];
        pts[3 * v + 1] = delaunay_->vertex_ptr(v)[1];
        pts[3 * v + 2] =
            (delaunay_->dimension() >= 3) ? delaunay_->vertex_ptr(v)[2] : 0.0;
    }

    if (delaunay_->keeps_infinite()) {

        // The convex hull can be efficiently traversed only if infinite
        // tetrahedra are kept.
        for (
            index_t t = delaunay_->nb_finite_cells();
            t < delaunay_->nb_cells(); ++t
            ) {
            for (index_t lv = 0; lv < 4; ++lv) {
                signed_index_t v = delaunay_->cell_vertex(t, lv);
                if (v != -1) {
                    tri2v.push_back(index_t(v));
                }
            }
        }

        mesh_.facets.assign_triangle_mesh(3, pts, tri2v, true);

        nbcell_ = delaunay_->nb_finite_cells();
    }
    else {

        nbcell_ = delaunay_->nb_cells();
    }

    if (delaunay_->dimension() == 3 + iweight_) {
        std::vector<index_t> tet2v(nbcell_ * 4);
        for (index_t t = 0; t < nbcell_; ++t) {
            tet2v[4 * t] = index_t(delaunay_->cell_vertex(t, 0));
            tet2v[4 * t + 1] = index_t(delaunay_->cell_vertex(t, 1));
            tet2v[4 * t + 2] = index_t(delaunay_->cell_vertex(t, 2));
            tet2v[4 * t + 3] = index_t(delaunay_->cell_vertex(t, 3));
        }
        mesh_.cells.assign_tet_mesh(3, pts, tet2v, true);
    } 
    else if (delaunay_->dimension() == 2 + iweight_) {
        tri2v.resize(nbcell_ * 3);
        for (index_t t = 0; t <nbcell_; ++t) {
            tri2v[3 * t] = index_t(delaunay_->cell_vertex(t, 0));
            tri2v[3 * t + 1] = index_t(delaunay_->cell_vertex(t, 1));
            tri2v[3 * t + 2] = index_t(delaunay_->cell_vertex(t, 2));
        }
        mesh_.facets.assign_triangle_mesh(3, pts, tri2v, true);
    }
    mesh_.show_stats();

    Logger::div("Saving the result");

    if (CmdLine::get_arg_bool("attributes")) {
        flags.set_attribute(MESH_FACET_REGION);
        flags.set_attribute(MESH_CELL_REGION);
    }
    if (GEO::CmdLine::get_arg_bool("single_precision")) {
        mesh_.vertices.set_single_precision();
    }
    if (FileSystem::extension(filename) == "geogram") {
        mesh_.vertices.set_double_precision();
    }
    bool result = true;
    if (mesh_save(mesh_, filename, flags)) {
        current_file_ = filename;
    }
    else {
        result = false;
    }

    return result;
}

The code and the output file are in the attachment. code+outfile.zip

[BrunoLevy]

Hello, In periodic mode, the same vertex can have several "instances", depending on where it is seen from. For instance, a vertex v1 near the side 'x=1' of the cube may be connected to another vertex v2 on the other side, that is, v2 is near the side x=0. When you consider v1 and its neighbors (for instance to compute the Voronoi cell), you need to translate v2 by [1 0 0]. This is encoded in the "instance" of v2: each vertex can be translated in 3x3x3 different ways. This explains why you see vertex indices that are (up to 27x) larger than the number of points. To get the (appropriately translated) geometry of a vertex, you need to use the PeriodicDelaunay3d::vertex() function instead of tri.vertex_ptr() that does not work for periodic vertices instances. If you do not need to translate the instance, you can use tri.vertex_ptr(tri.periodic_vertex_real(v)) instead (where tri is your PeriodicDelaunayTriangulation).

[luogyong]

Hi Bruno, thank you for your prompt reply.

Your recommended approach, Periodically3d::vertex(), helped me obtain the correct vertex coordinates. However, I am unsure how to ensure that the vertex indices of a Tet element fall within the range of 1 to Nb_vertices (410). The total number of vertices output is around 410 (with 100 non-periodic points), yet the Tet vertex indices generated range from 1 to 2700.

Thank you.

[BrunoLevy]

Use PeriodicDelaunayTriangulation3d::periodic_vertex_real(v) to transform a vertex instance into the associated real vertex.

[luogyong]

Finally, I resolved this problem by implementing a linking array that maps a periodic vertex to the corresponding output nodal index. The save code is as shown below:

    bool DemoDelaunay3dApplication::save(const std::string& filename) {
        MeshIOFlags flags;

        vector<double> pts;
        vector<index_t> tet2v,pv2outputindices(delaunay_->nb_vertices_non_periodic_*27);

        if (delaunay_->keeps_infinite()) {

            nbcell_ = delaunay_->nb_finite_cells();
        }
        else {

            nbcell_ = delaunay_->nb_cells();
        }

        index_t cellsize;

        if (delaunay_->dimension() == 3 + iweight_) {
            cellsize = 4;

            //mesh_.cells.assign_tet_mesh(3, pts, tet2v, true);

        }
        else if (delaunay_->dimension() == 2 + iweight_) {
            cellsize = 3;
            //mesh_.facets.assign_triangle_mesh(3, pts, tri2v, true);
        }

        tet2v.resize(nbcell_ * cellsize, 0);
        for (index_t t = 0; t < nbcell_; ++t) {
            for (index_t i = 0; i < cellsize; ++i) {
                tet2v[cellsize * t + i] = index_t(delaunay_->cell_vertex(t, i));
                pv2outputindices[tet2v[cellsize * t + i]] = 1;
            }
        }

        index_t n = count_if(pv2outputindices.begin(), pv2outputindices.end(), [](index_t t) {return t == 1; });
        pts.resize(3*n);
        index_t j = 0;
        for (index_t i = 0; i < pv2outputindices.size(); ++i) {
            if (pv2outputindices[i] > 0) {
                pts[3 * j] = delaunay_->vertex(i)[0];
                pts[3 * j+1] = delaunay_->vertex(i)[1];
                pts[3 * j+2] = (delaunay_->dimension() >= 3) ? delaunay_->vertex(i)[2] : 0.0;
                pv2outputindices[i] = j++;
            }
        }

        for (index_t i = 0; i < tet2v.size(); ++i) {
            tet2v[i]=pv2outputindices[tet2v[i]];
        }

        if (delaunay_->dimension() == 3 + iweight_) {

            mesh_.cells.assign_tet_mesh(3, pts, tet2v, true);

        }
        else if (delaunay_->dimension() == 2 + iweight_) {

            mesh_.facets.assign_triangle_mesh(3, pts, tet2v, true);
        }

        if (delaunay_->keeps_infinite()) {

            // The convex hull can be efficiently traversed only if infinite
        // tetrahedra are kept.
            //geo_assert(delaunay_->keeps_infinite());

            // The convex hull can be retrieved as the finite facets
            // of the infinite cells (note: it would be also possible to
            // throw away the infinite cells and get the convex hull as
            // the facets adjacent to no cell). Here we use the infinite
            // cells to show an example with them.

            //// This block is just a sanity check
            //{
            //  for (index_t t = 0; t < delaunay_->nb_finite_cells(); ++t) {
            //      geo_debug_assert(delaunay_->cell_is_finite(t));
            //  }

            //  for (index_t t = delaunay_->nb_finite_cells();
            //      t < delaunay_->nb_cells(); ++t) {
            //      geo_debug_assert(delaunay_->cell_is_infinite(t));
            //  }
            //}

            // This iterates on the infinite cells
            tet2v.clear();
            for (
                index_t t = delaunay_->nb_finite_cells();
                t < delaunay_->nb_cells(); ++t
                ) {
                for (index_t lv = 0; lv < 4; ++lv) {
                    signed_index_t v = delaunay_->cell_vertex(t, lv);
                    if (v != -1) {
                        tet2v.push_back(pv2outputindices[index_t(v)]);
                    }
                }
            }

            mesh_.facets.assign_triangle_mesh(3, pts, tet2v, true);

            //nbcell_ = delaunay_->nb_finite_cells();
        }

        mesh_.show_stats();

        Logger::div("Saving the result");

        if (CmdLine::get_arg_bool("attributes")) {
            flags.set_attribute(MESH_FACET_REGION);
            flags.set_attribute(MESH_CELL_REGION);
        }
        if (GEO::CmdLine::get_arg_bool("single_precision")) {
            mesh_.vertices.set_single_precision();
        }
        if (FileSystem::extension(filename) == "geogram") {
            mesh_.vertices.set_double_precision();
        }
        bool result = true;
        if (mesh_save(mesh_, filename, flags)) {
            current_file_ = filename;
        }
        else {
            result = false;
        }

        return result;
    }