小贝子编程

如何防止CGAL在三维Alpha形状中生成简并四面体



我正试图在三维中检索阿尔法形状内的四面体，以进行一些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的输出为(以gmsh表示(：

CGAL在gmsh 中的输出

然而，图中似乎有一些交叉/重叠的元素，当程序在调试中通过gdb执行时，在alpha形状构造函数处有一个断点，会抛出以下内容：文件Undefine.h:中的"CGAL:：Undefinable conversion">

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;
}

如果我稍微扰动点的位置(在我的程序中，我使用gmsh生成初始状态，节点位置不是完全"完整"的(，就会出现那些奇怪的元素，而如果节点位置是完全整数，则元素的数量是正确的(在这种情况下为1017(。

好吧，我找到了一种防止这些碎片出现的方法。我在gmsh中生成"；超限；网格在边界上，这意味着节点的位置就像在一个完美的正方形网格中一样。通过放松这个条件，我最终得到了一组在边界上形成等边三角形的初始节点，这些节点不会生成切片。

我仍然需要找到一种方法来防止网格变形时出现碎片，但这与CGAL无关。

您的代码基本上由位于规则网格上的点组成，然后您会稍微扰动这些点。结果，你最终得到了通常被称为slivers的东西(在你最喜欢的搜索引擎中寻找"slivers tet mesh"(。我在这里的建议是使用规则网格生成初始网格，并在输入上应用扰动(由于顶点上有id，因此可以很容易地将扰动存储在循环中(。

相关内容

最新更新


热门标签：

javascript python java c# php android html jquery c++ css ios sql mysql arrays asp.net json python-3.x ruby-on-rails .net sql-server django objective-c excel regex ruby linux ajax iphone xml vba spring asp.net-mvc database wordpress string postgresql wpf windows xcode bash git oracle list vb.net multithreading eclipse algorithm macos powershell visual-studio image forms numpy scala function api selenium