search

CGAL / cgal

The public CGAL repository, see the README below
https://github.com/CGAL/cgal#readme
Other
4.87k stars 1.38k forks source link

representation of hexahedral mesh using linear cell complex #6364

Open Supranaturaler opened 2 years ago

Supranaturaler commented 2 years ago

Issue Details

Hi guys, I want to use linear cell complex to represent hexahedral mesh. Below is the mesh I want to construct. image I have all the coordinates in each vertex and I'm using Linear_cell_complex_incremental_builder_3 to construct LCC. When I run this code in debug mode, CGAL assertion fails. How can I build the mesh in a proper way?

CGAL error: assertion violation! Expression : lcc.template is_free<2>(opposite) File : \3rdparty\CGAL-5.4\include\CGAL\Linear_cell_complex_incremental_builder.h Line : 119 Explanation: Refer to the bug-reporting instructions at https://www.cgal.org/bug_report.html Describe your issue. Please be specific (compilation error, runtime error, unexpected behavior, wrong results, etc.).

Source Code

#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;

#include "CGAL/Linear_cell_complex_for_combinatorial_map.h"
typedef CGAL::Linear_cell_complex_traits<3, Kernel> LCCTraits;
typedef CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, LCCTraits> LCCMesh;

#include "CGAL/Linear_cell_complex_constructors.h"
typedef CGAL::Linear_cell_complex_incremental_builder_3<LCCMesh> LCCBuilder;

int main(int argc, char* argv[])
{
    std::vector<double> vecX = { 0,1,2,3,4,5,6,7,8 };
    int xSize = vecX.size();
    std::vector<double> vecY = { 0,1,2,3,4,5,6,7,8 };
    int ySize = vecY.size();
    std::vector<double> vecZ = { 0,1,2,3,4,5,6,7,8 };
    int zSize = vecZ.size();

    LCCMesh lccMesh;
    LCCBuilder lccBuilder(lccMesh);
    int facetSize = xSize * (ySize - 1) * (zSize - 1) +
        ySize * (xSize - 1) * (zSize - 1) + zSize * (xSize - 1) * (ySize - 1);
    lccBuilder.begin_surface(xSize * ySize * zSize, facetSize, 0);
    for (int i = 0; i < xSize; ++i)
    {
        for (int j = 0; j < ySize; ++j)
        {
            for (int k = 0; k < zSize; ++k)
            {
                lccBuilder.add_vertex(Point_3(vecX[i], vecY[j], vecZ[k]));
            }
        }
    }

    int startIndex = 0;
    for (int i = 0; i < xSize; ++i)
    {
        for (int j = 0; j < ySize; ++j)
        {
            for (int k = 0; k < zSize; ++k)
            {
                if (j != ySize - 1 && k != zSize - 1)
                {
                    int index1 = startIndex;
                    int index2 = startIndex + zSize;
                    int index3 = startIndex + zSize + 1;
                    int index4 = startIndex + 1;

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index1);
                    lccBuilder.add_vertex_to_facet(index2);
                    lccBuilder.add_vertex_to_facet(index3);
                    lccBuilder.add_vertex_to_facet(index4);
                    lccBuilder.end_facet();
                }

                startIndex++;
            }
        }
    }

    startIndex = 0;
    for (int i = 0; i < xSize; ++i)
    {
        for (int j = 0; j < ySize; ++j)
        {
            for (int k = 0; k < zSize; ++k)
            {
                if (i !=xSize - 1 && k != zSize - 1)
                {
                    int index1 = startIndex;
                    int index2 = startIndex + 1;
                    int index3 = startIndex + ySize * zSize + 1;
                    int index4 = startIndex + ySize * zSize;

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index1);
                    lccBuilder.add_vertex_to_facet(index2);
                    lccBuilder.add_vertex_to_facet(index3);
                    lccBuilder.add_vertex_to_facet(index4);
                    lccBuilder.end_facet();
                }

                startIndex++;
            }
        }
    }

    startIndex = 0;
    for (int i = 0; i < xSize; ++i)
    {
        for (int j = 0; j < ySize; ++j)
        {
            for (int k = 0; k < zSize; ++k)
            {
                if (i != xSize - 1 && j != ySize - 1)
                {
                    int index1 = startIndex;
                    int index2 = startIndex + ySize * zSize;
                    int index3 = startIndex + ySize * zSize + zSize;
                    int index4 = startIndex + zSize;

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index1);
                    lccBuilder.add_vertex_to_facet(index2);
                    lccBuilder.add_vertex_to_facet(index3);
                    lccBuilder.add_vertex_to_facet(index4);
                    lccBuilder.end_facet();
                }

                startIndex++;
            }
        }
    }

    lccBuilder.end_surface();
    return EXIT_SUCCESS;
}

Environment

gdamiand commented 2 years ago

1) You need to do a lccBuilder.start_surface() and lccBuilder.end_surface(); for each cube. Indeed, each cube is described by one different surface.

2) Orientation must be globally consistent: the 4 vertices of the face of one cube must be in reverse order for the same face of the other cube incident to the same face.

gdamiand commented 2 years ago

One more remark: you need to use cgal master because the functionality allowing to use the incremental builder to create different volumes was introduced in #6055 (and thus is not in cgal 5.4).

Supranaturaler commented 2 years ago
  1. You need to do a lccBuilder.start_surface() and lccBuilder.end_surface(); for each cube. Indeed, each cube is described by one different surface.

I want the adjacent two cubes to share one face. So there will be 9 9 9 vertices in total. If I use lccBuilder.start_surface() and lccBuilder.end_surface() for each cube, won't there be 8 8 8 * 8 vertices? Many of them are duplicate.

gdamiand commented 2 years ago

I want the adjacent two cubes to share one face. So there will be 9 9 9 vertices in total. If I use lccBuilder.start_surface() and lccBuilder.end_surface() for each cube, won't there be 8 8 8 * 8 vertices? Many of them are duplicate.

The vertices will not be duplicated. You start to call add_vertex before any start_surface. All these vertices will be shared between all the surfaces.

Supranaturaler commented 2 years ago

The function begin_surface() has been modified in #6055. It will be more convenient to construct LCC. Thank you for your reply! I will try this improvement in master branch.

Supranaturaler commented 2 years ago

@gdamiand Hi, I test the new code in master branch. There seems to be a bug in CGAL::write_off(LCC& alcc, std::ostream& out). The face size is not correct. Please check.

#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;

#include "CGAL/Linear_cell_complex_for_combinatorial_map.h"
typedef CGAL::Linear_cell_complex_traits<3, Kernel> LCCTraits;
typedef CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, LCCTraits> LCCMesh;

#include "CGAL/Linear_cell_complex_constructors.h"
typedef CGAL::Linear_cell_complex_incremental_builder_3<LCCMesh> LCCBuilder;

#include "CGAL/draw_linear_cell_complex.h"

int main(int argc, char* argv[])
{
    std::vector<double> vecX = { 0,1,2,3,4,5,6,7,8 };
    int xSize = vecX.size();
    std::vector<double> vecY = { 0,1,2,3,4,5,6,7,8 };
    int ySize = vecY.size();
    std::vector<double> vecZ = { 0,1,2,3,4,5,6,7,8 };
    int zSize = vecZ.size();

    LCCMesh lccMesh;
    LCCBuilder lccBuilder(lccMesh);

    for (int i = 0; i < xSize; ++i)
    {
        for (int j = 0; j < ySize; ++j)
        {
            for (int k = 0; k < zSize; ++k)
            {
                lccBuilder.add_vertex(Point_3(vecX[i], vecY[j], vecZ[k]));
            }
        }
    }

    int startIndex = 0;
    for (int i = 0; i < xSize; ++i)
    {
        for (int j = 0; j < ySize; ++j)
        {
            for (int k = 0; k < zSize; ++k)
            {
                if (i != xSize - 1 && j != ySize - 1 && k != zSize - 1)
                {     
                    lccBuilder.begin_surface();
                    int index1 = startIndex;
                    int index2 = startIndex + zSize;
                    int index3 = startIndex + ySize * zSize + zSize;
                    int index4 = startIndex + ySize * zSize;
                    int index5 = startIndex + ySize * zSize + 1;
                    int index6 = startIndex + 1;
                    int index7 = startIndex + zSize + 1;
                    int index8 = startIndex + ySize * zSize + zSize + 1;

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index1);
                    lccBuilder.add_vertex_to_facet(index2);
                    lccBuilder.add_vertex_to_facet(index3);
                    lccBuilder.add_vertex_to_facet(index4);
                    lccBuilder.end_facet();

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index5);
                    lccBuilder.add_vertex_to_facet(index8);
                    lccBuilder.add_vertex_to_facet(index7);
                    lccBuilder.add_vertex_to_facet(index6);
                    lccBuilder.end_facet();

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index1);
                    lccBuilder.add_vertex_to_facet(index4);
                    lccBuilder.add_vertex_to_facet(index5);
                    lccBuilder.add_vertex_to_facet(index6);
                    lccBuilder.end_facet();

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index3);
                    lccBuilder.add_vertex_to_facet(index2);
                    lccBuilder.add_vertex_to_facet(index7);
                    lccBuilder.add_vertex_to_facet(index8);
                    lccBuilder.end_facet();

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index5);
                    lccBuilder.add_vertex_to_facet(index4);
                    lccBuilder.add_vertex_to_facet(index3);
                    lccBuilder.add_vertex_to_facet(index8);
                    lccBuilder.end_facet();

                    lccBuilder.begin_facet();
                    lccBuilder.add_vertex_to_facet(index1);
                    lccBuilder.add_vertex_to_facet(index6);
                    lccBuilder.add_vertex_to_facet(index7);
                    lccBuilder.add_vertex_to_facet(index2);
                    lccBuilder.end_facet();

                    lccBuilder.end_surface();
                }

                startIndex++;
            }
        }
    }

    CGAL::draw(lccMesh);
    std::ofstream stream;
    stream.open("output.off");
    CGAL::write_off(lccMesh, stream);
    stream.close();

    return EXIT_SUCCESS;
}

There should be 8 8 8 * 6 = 3072 faces. But in the export off file, the face size is 1728. image image

gdamiand commented 2 years ago

off file format is not able to represent volumic meshes. And your computation of number of faces is wrong: between two adjacent volumes, there is only 1 face (shared by the two cubes) and not 2.

Supranaturaler commented 2 years ago

If so, the function CGAL::write_off(LCC& alcc, std::ostream& out) should be modified. The face size is not match with its list of faces. Face size is 1728, but the list of faces is 3072.

image image