Как предотвратить генерирование CGAL вырожденного тетраэдра в 3D Alpha Shape? - PullRequest
Новая
       12

Как предотвратить генерирование CGAL вырожденного тетраэдра в 3D Alpha Shape?

0 голосов
/ 30 марта 2020

Я пытаюсь получить тетраэдры внутри альфа-формы в 3D для некоторых вычислений FEM. Код, который я использую, выглядит следующим образом: 

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>


typedef CGAL::Exact_predicates_inexact_constructions_kernel                 Kernel;
typedef Kernel::FT                                                          FT;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::size_t, Kernel>    Vb3;
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb3>                        asVb3;
typedef CGAL::Alpha_shape_cell_base_3<Kernel>                               asCb3;
typedef CGAL::Triangulation_data_structure_3<asVb3,asCb3>                   asTds3;
typedef CGAL::Delaunay_triangulation_3<Kernel,asTds3, CGAL::Fast_location>  asTriangulation_3;
typedef CGAL::Alpha_shape_3<asTriangulation_3>                              Alpha_shape_3;
typedef Kernel::Point_3                                                     Point_3;
typedef Kernel::Tetrahedron_3                                               Tetrahedron_3;

void triangulateAlphaShape3D()
{
    ...

    const Alpha_shape_3 as(pointsList.begin(), pointsList.end(),
                           FT(alpha),
                           Alpha_shape_3::GENERAL);

    for(Alpha_shape_3::Finite_cells_iterator fit = as.finite_cells_begin() ;
        fit != as.finite_cells_end() ; ++fit)
    {
        if(as.classify(fit) == Alpha_shape_3::INTERIOR)
        {
            const Alpha_shape_3::Cell_handle cell{fit};

            const std::vector<std::size_t> element{cell->vertex(0)->info(),
                                                   cell->vertex(1)->info(),
                                                   cell->vertex(2)->info(),
                                                   cell->vertex(3)->info()};
        }
    }

    ...
}

Мой список точек состоит из точек, которые образуют куб в блоке. Вывод CGAL (визуализирован в gm sh):

Вывод CGAL визуализирован в gm sh

Однако кажется, что есть какое-то пересечение / перекрывающиеся элементы на рисунке, и программа выдает следующее при выполнении через gdb при отладке с точкой останова в конструкторе формы альфа: «Неразрешимое преобразование CGAL :: Untermin» из файла Untermin.h: 

T make_certain() const
  {
    if (is_certain())
      return _i;
    CGAL_PROFILER(std::string("Uncertain_conversion_exception thrown for CGAL::Uncertain< ")
          + typeid(T).name() + " >");
    throw Uncertain_conversion_exception(
                  "Undecidable conversion of CGAL::Uncertain<T>");
  }

Стек вызовов: 

#0 0xa31c08 libstdc++-6!.cxa_throw() (libstdc++-6.dll:??)
#1 0x62060a4a   CGAL::Uncertain<bool>::make_certain(this=0x22d0ac) (CGAL/Uncertain.h:125)
#2 0x62060ae5   CGAL::Uncertain<bool>::operator bool(this=0x22d0ac) (CGAL/Uncertain.h:133)
#3 0x620240bd   CGAL::coplanar_orientationC3<CGAL::Interval_nt<false> >(px=..., py=..., pz=..., qx=..., qy=..., qz=..., rx=..., ry=..., rz=...) (CGAL/predicates/kernel_ftC3.h:209)
#4 0x62056e68   CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >::operator() (this=0x22d5af, p=..., q=..., r=...) (CGAL/Cartesian/function_objects.h:3649)
#5 0x620518ea   CGAL::Filtered_predicate<CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Mpzf> >, CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >, CGAL::Cartesian_converter<CGAL::Type_equality_wrapper<CGAL::Cartesian_base_no_ref_count<double, CGAL::Epick>, CGAL::Epick>, CGAL::Simple_cartesian<CGAL::Mpzf>, CGAL::NT_converter<double, CGAL::Mpzf> >, CGAL::Cartesian_converter<CGAL::Type_equality_wrapper<CGAL::Cartesian_base_no_ref_cou (CGAL/Filtered_predicate.h:167)
#6 0x6205081f   CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:589)
#7 0x62059ec2   CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1509)
#8 0x6205950e   CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1545)
#9 0x6205988c   CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:775)
#10 0x620599c9  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:857)
#11 0x6204e578  CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:1323)
#12 0x6201d924  CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:3696)
#13 0x620254cc  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1165)
#14 0x620257bd  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1150)
#15 0x62025636  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:537)
#16 0x62027f9c  CGAL::Triangulation_hierarchy_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CG (CGAL/Triangulation_hierarchy_3.h:279)
#17 0x6202816f  CGAL::Triangulation_hierarchy_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CG (CGAL/Triangulation_hierarchy_3.h:300)
#18 0x62017b27  CGAL::Alpha_shape_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Sequential_tag>, CGAL (CGAL/Alpha_shape_3.h:269)

Это какие-то вырожденные тетраэдры? Как я могу предотвратить генерацию их CGAL?

В настоящее время я использую CGAL 4.11.

Обновление : я перешел на CGAL 5.0.2 и проблема с броском исчезла. Однако у меня все еще есть эта странная проблема. Тестовый код: 

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>

#include <gmsh.h>

#include <algorithm>
#include <fstream>
#include <random>

typedef CGAL::Exact_predicates_inexact_constructions_kernel                 Kernel;
typedef Kernel::FT                                                          FT;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::size_t, Kernel>    Vb3;
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb3>                        asVb3;
typedef CGAL::Alpha_shape_cell_base_3<Kernel>                               asCb3;
typedef CGAL::Triangulation_data_structure_3<asVb3,asCb3>                   asTds3;
typedef CGAL::Delaunay_triangulation_3<Kernel,asTds3, CGAL::Fast_location>  asTriangulation_3;
typedef CGAL::Alpha_shape_3<asTriangulation_3>                              Alpha_shape_3;
typedef Kernel::Point_3                                                     Point_3;
typedef Kernel::Tetrahedron_3                                               Tetrahedron_3;

int main(int argc, char **argv)
{
    std::random_device rd;
    std::mt19937 mt(rd());
    std::uniform_real_distribution<double> dist(-0.000000001, 0.000000001);
    //std::uniform_real_distribution<double> dist(0,0);

    std::vector<std::vector<std::size_t>> elementsList;

    const double alpha = 0.9;

    // We have to construct an intermediate representation for CGAL. We also reset
    // nodes properties.
    std::vector<std::pair<Point_3, std::size_t>> pointsList;

    std::size_t counter = 0;

    /********************************************************************************
                                     Box Nodes
     *******************************************************************************/
    for(double y = 1; y < 10 ; y += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(10, y + dist(mt), z + dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double x = 0; x <= 10 ; x += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), 0, z+ dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double y = 1; y <= 10 ; y += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(0, y+ dist(mt), z+ dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double x = 1; x <= 10 ; x += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), 10, z+ dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double x = 0; x <= 10 ; x += 1)
    {
        for(double y = 0; y <= 10 ; y += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), y+ dist(mt), 0), counter);
            pointsList.push_back(point);
            counter++;
        }
    }


    /********************************************************************************
                                     Cube nodes
     *******************************************************************************/
    for(double x = 1; x <= 5 ; x += 1)
    {
        for(double y = 1; y <= 5 ; y += 1)
        {
            for(double z = 1; z <= 5 ; z += 1)
            {
                std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), y+ dist(mt), z+ dist(mt)), counter);
                pointsList.push_back(point);
                counter++;
            }
        }
    }

    std::random_shuffle(pointsList.begin(), pointsList.end());

    for(std::size_t n = 0; n < pointsList.size() ; ++n)
    {
        pointsList[n].second = n;
    }

    for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
    {
        for(std::size_t n2 = 0 ; n2 < pointsList.size() ; ++n2)
        {
            if(n != n2)
            {
                if(pointsList[n].first == pointsList[n2].first)
                {
                    std::cerr << "Duplicate node found: " << "(" << pointsList[n].first.x() << ", "
                                                                 << pointsList[n].first.y() << ", "
                                                                 << pointsList[n].first.z() << ")";
                    std::cerr << std::endl;
                    std::cerr << "Counter: " << pointsList[n].second << " vs " << pointsList[n2].second << std::endl;
                    throw std::runtime_error("Duplicates nodes!");
                }
            }
        }
    }

    std::cout << "Number of point: " << pointsList.size() << std::endl;

    const Alpha_shape_3 as(pointsList.begin(), pointsList.end(),
                           FT(alpha),
                           Alpha_shape_3::GENERAL);

    // We check for each triangle which one will be kept (alpha shape), then we
    // perfom operations on the remaining elements
    for(Alpha_shape_3::Finite_cells_iterator fit = as.finite_cells_begin() ;
        fit != as.finite_cells_end() ; ++fit)
    {
        // If true, the elements are fluid elements
        if(as.classify(fit) == Alpha_shape_3::INTERIOR)
        {
            const Alpha_shape_3::Cell_handle cell{fit};

            const std::vector<std::size_t> element{cell->vertex(0)->info(),
                                                   cell->vertex(1)->info(),
                                                   cell->vertex(2)->info(),
                                                   cell->vertex(3)->info()};

            elementsList.push_back(element);
        }
    }

    std::cout << "Number of Element: " << elementsList.size() << std::endl;

    gmsh::initialize();
    gmsh::model::add("theModel");
    gmsh::model::setCurrent("theModel");
    gmsh::model::addDiscreteEntity(0, 1);
    gmsh::model::addDiscreteEntity(3, 2);
    gmsh::view::add("Nothing", 1);

    std::vector<std::size_t> nodesTags(pointsList.size());
    std::vector<double> nodesCoord(3*pointsList.size());

    for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
    {
        nodesTags[n] = n + 1;
        nodesCoord[3*n] = pointsList[n].first.x();
        nodesCoord[3*n + 1] = pointsList[n].first.y();
        nodesCoord[3*n + 2] = pointsList[n].first.z();
    }

    gmsh::model::mesh::addNodes(0, 1, nodesTags, nodesCoord);

    std::vector<std::size_t> elementTags(elementsList.size());
    std::vector<std::size_t> nodesTagsPerElement(4*elementsList.size());
    for(std::size_t elm = 0 ; elm < elementsList.size() ; ++elm)
    {
        elementTags[elm] = elm + 1;
        for(std::size_t n = 0 ; n < 4 ; ++n)
            nodesTagsPerElement[4*elm + n] = pointsList[elementsList[elm][n]].second + 1;
    }

    gmsh::model::mesh::addElementsByType(2, 4, elementTags, nodesTagsPerElement); //Tetrahedron 4
    gmsh::model::mesh::addElementsByType(1, 15, nodesTags, nodesTags); //Point

    std::vector<std::vector<double>> dataNothing(pointsList.size());

    for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
    {
        const std::vector<double> nothingVec{0};
        dataNothing[n] = nothingVec;
    }

    gmsh::view::addModelData(1, 0, "theModel", "NodeData", nodesTags, dataNothing, 0, 1);

    gmsh::view::write(1, "mesh.msh", false);

    gmsh::model::remove();

    gmsh::finalize();

    return 0;
}

Если я слегка возмущаю положение точек (в моей программе я генерирую исходное состояние, используя gm sh, а положение узла не совсем "целое"), то этот странный элемент появляется, пока если позиции узлов являются целыми числами, номер элемента правильный (в данном случае 1017).

1 Ответ

0 голосов
/ 02 апреля 2020

Ваш код в основном состоит из точек, находящихся на регулярной сетке, которые вы затем слегка возмущаете В результате вы получаете то, что обычно называют «мычками» (ищите «любезно меня» sh в вашей любимой поисковой системе). Мой совет здесь заключается в том, чтобы сгенерировать начальное значение me sh с использованием регулярной сетки и применить возмущение к входу (вы можете легко сохранить возмущение в вашем l oop, поскольку у вас есть идентификаторы на вершинах).

...