希望基于CGAL，写一个自己的网格生成与优化程序。 先从2D的情形开始吧！
     首先从功能谈起。 这个程序要以文件形式输入区域信息，能够处理曲边的情形， 生成满足一定密度的Delaunay网格，作为一个初始的的有限元计算网格；然后，给定一个密度函数（连续的或者离散的）， 再调用CVT，ODT方法优化网格，不断的重复这个过程序，直到满足一定的条件为止。

* 区域描述文件的数据格式：

采用sample.poly的数据格式，如下：

.poly files

• First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>
• Following lines: <vertex #> <x> <y> [attributes] [boundary marker]
• One line: <# of segments> <# of boundary markers (0 or 1)>
• Following lines: <segment #> <endpoint> <endpoint> [boundary marker]
• One line: <# of holes>
• Following lines: <hole #> <x> <y>
• One line: <# of regional attributes 
• Following lines: <region #> <x> <y>
  

  # square.poly
   4 2 0 1
     1   -1 -1   1
     2   1 -1   1
     3   1 1   1
     4   -1 1   1
     4  1
     1   1 2   3
     2   2 3   4
     3   3 4   5
     4   4 1   6
     0

#2hexa_in_square.poly
     16     2     0     1
     1        0.0000000000        0.0000000000     1
     2        1.0000000000        0.0000000000     1
     3        1.0000000000        1.0000000000     1
     4        0.0000000000        1.0000000000     1
     5        0.1500000060        0.7500000000     1
     6        0.2002999932        0.6633999944     1
     7        0.3001999855        0.6639999747     1
     8        0.3499999940        0.7500000000     1
     9        0.3001999855        0.8360000253     1
    10        0.2002999932        0.8366000056     1
    11        0.5000000000        0.4000000060     1
    12        0.5503000021        0.3134000003     1
    13        0.6502000093        0.3140000105     1
    14        0.6999999881        0.4000000060     1
    15        0.6502000093        0.4860000014     1
    16        0.5503000021        0.4866000116     1
    16 1
    1    1 2    3
    2    2 3    4   
    3    3 4    5
    4     4 1    6
    5    5 6    7
    6    6 7    8
    7    7 8    9
    8    8 9    10
    9    9 10    11
    10    10 5    12
    11    11 12    13
    12    12 13    14
    13    13 14    15
    14    14 15    16
    15    15 16    17
    16    16 11    18
     2
     1        0.2500000000        0.7500000000
     2        0.6000000238        0.4000000060

如果想处理拥有曲边的区域情形， 首先要有一个对区域的多边形逼近， 再加上一个距离函数，能够把新生成的边界点投影到边界上。

在.poly文件当中，边界点的标记必须是1或2， 如果是1，表示它在整个网格生成过程中是固定的；如果是2，则在直边的情形下它是固定的，在曲边的情形下，它要被投回到边界上。

顶点和面要加入额外的信息。顶点要知道它是否在边界上，是否要投影回边界，对应一个尺寸函数值。面要知道它属于那个子区域。

* 算法的实现CGAL基础
    ＊读入区域文件。
    ＊设定限止边和种子
    ＊进行初始的剖分，形成一个均匀的网格
    ＊如果有曲边，把本就是曲边的点投影到曲边上。

     ＊根据密度函数插入新的节点，形成一个满足一定密度函数的网格

     ＊网格优化：CVT， ODT， distmesh

2， 输出数据的文件格式

3，可视化接口的实现

利用Geoview 来实现可视化。

2010.09.13

 昨天网格生成的程序可以运行了，在此基础之上就能进行更深入的学习和实验。

在用std::vector时，我先给这种类型的变量Vhs的长度设为n，然后又用push_back把读入的数据存在变量当中（实际上是存在n＋1的位置上）。但我以为是存开头，用Vhs[0]从而导致访问的错误。最终的程序如下：


#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesher_2.h>

#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Constrained_triangulation_face_base_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>

#include <CGAL/Triangulation_vertex_base_with_info_2.h>



#include <iostream>
#include <fstream>

#include <unistd.h> // for sleep()


#include <CGAL/IO/Geomview_stream.h>
#include <CGAL/IO/Triangulation_geomview_ostream_2.h>

#include <CGAL/Bbox_2.h>
#include <CGAL/Polygon_2_algorithms.h>

class Face_info
{
  int attribute;
public:
  Face_info():attribute(1){}
  Face_info(int a): attribute(a) {}
  int get_attribute()
  {
    return attribute;
  }
  void set_attribute(int a)
  {
    attribute = a;
  }

};

class Vertex_info
{
  int bd_flag;// 1:fixed point; 2:should project into boundary; 0: other points. 
public:
  Vertex_info():bd_flag(0){}
  Vertex_info(int bf):bd_flag(bf){}
  void set_bd_flag(int bf)
  {
    bd_flag = bf;
  }
  int get_bd_flag()
  {
    return bd_flag;
  }
};

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Triangulation_vertex_base_with_info_2<Vertex_info,K> Vb;


typedef CGAL::Triangulation_face_base_with_info_2<Face_info,K> Fbbb;
typedef CGAL::Constrained_triangulation_face_base_2<K,Fbbb> Fbb;
typedef CGAL::Delaunay_mesh_face_base_2<K,Fbb> Fb;


typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
typedef CGAL::Delaunay_mesher_2<CDT, Criteria> Mesher;

typedef CDT::Vertex_handle Vertex_handle;
typedef CDT::Face_handle Face_handle;
typedef CDT::Point Point;


typedef std::pair<int, int>   Myedge;
typedef std::list<Myedge> List_edges;
typedef std::vector<Point> Vector_points;

int main()
{
  CDT cdt;

  // begin of reading .poly file

  int n_point=0;
  Vector_points con_points;
  std::vector<int> bd_flag;

  int n_con_edge=0;
  List_edges con_edges;

  int n_seeds;
  std::list<Point> seeds;

  int n_attri=0;
  Vector_points attri_points;
  std::vector<int> attri_flag;

  double x,y;
  Myedge e;
  int temp;

  std::cout<<"begin to read!"<<std::endl;

  std::ifstream input("squa_in_squa.poly",std::ios::in);

  input>>n_point>>temp>>temp>>temp;

  std::cout<<"n_point="<<n_point<<std::endl;


  for(int i=0;i<n_point;i++)
    {
      input>>temp>>x>>y>>temp;
      con_points.push_back(Point(x,y));
      std::cout<<x<<","<<y<<std::endl;
      bd_flag.push_back(temp);
    }

  std::cout<<"finished reading bounadary points!"<<std::endl;
 
  input>>n_con_edge>>temp;

 

  for(int i=0;i<n_con_edge;i++)
    {
      input>>temp>>e.first>>e.second>>temp;
      con_edges.push_back(e);
    }

  std::cout<<"n_con_edge="<<n_con_edge<<std::endl;
  input>>n_seeds;

  std::cout<<"n_seeds="<<n_seeds<<std::endl;

  if(n_seeds)
    {
      for(int i=0;i<n_seeds;i++)
    {
      input>>temp>>x>>y;
      seeds.push_back(Point(x,y));
    }
      input>>n_attri;
      if(n_attri)
    {
      for(int i=0;i<n_attri;i++)
        {
          input>>temp>>x>>y>>temp;
          attri_points.push_back(Point(x,y));
          attri_flag.push_back(temp);
        }
    }
    }
  input.close();

  // end of reading

  std::cout<<"computing the 2d box including the domain"<<std::endl;

  //computing the 2d box including the domain

  Vector_points::iterator left = CGAL::left_vertex_2(con_points.begin(),con_points.end(), K());
  Vector_points::iterator right = CGAL::right_vertex_2(con_points.begin(),con_points.end(),K());
  Vector_points::iterator top = CGAL::top_vertex_2(con_points.begin(),con_points.end(),K());
  Vector_points::iterator bottom = CGAL::bottom_vertex_2(con_points.begin(),con_points.end(),K());
  double w = (left->x() - right->x())/50.0;
  double h = (top->y() - bottom->y())/50.0;




  //the end of computing

  std::vector<Vertex_handle> Vhs;

 
  for(int i=0;i<n_point;i++)
    {
      Vertex_handle va=cdt.insert(con_points[i]);
      Vhs.push_back(va);

      std::cout<<"why!"<<"i="<<i<<std::endl;

      Vhs[i]->info().set_bd_flag(bd_flag[i]);
    }
 
  List_edges::iterator it;
  for(it = con_edges.begin();it!=con_edges.end();it++)
    {
      cdt.insert_constraint(Vhs[it->first-1],Vhs[it->second-1]);
    }


  std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl;

  std::cout << "Meshing the triangulation with default criterias..."
            << std::endl;

  Mesher mesher(cdt);
  if(n_seeds)
    {
      mesher.set_seeds(seeds.begin(),seeds.end());
      std::cout<<n_seeds<<std::endl;
    }

  mesher.set_criteria(Criteria(0.125, 0.5));
  mesher.refine_mesh();

  std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl;


   // Visualization

  CGAL::Geomview_stream gv(CGAL::Bbox_3(left->x()-w, bottom->y()-h, 0, right->x()+w, top->y()+h, 0));
  gv.set_line_width(4);
  gv.set_bg_color(CGAL::Color(0, 200, 200));

  std::cout << "Drawing 2D Delaunay triangulation in non-wired mode.n";
  gv.set_wired(true);
  gv << cdt;
  sleep(10);

  std::cout << "Enter a key to finish" << std::endl;
  char ch;
  std::cin >> ch;
  return 0;
}


转载本文请联系原作者获取授权，同时请注明本文来自魏华祎科学网博客。
链接地址：https://blog.sciencenet.cn/blog-284809-361796.html

上一篇：留学准备
下一篇：我的Ubuntu 9.10 软件配置