← Back to branch summary

~corrado-maurini/dolfin/tao

« back to all changes in this revision

Viewing changes to dolfin/generation/CSGCGALMeshGenerator2D.cpp

  • Committer: corrado maurini
  • Date: 2012-12-18 12:16:08 UTC
  • mfrom: (6685.78.207 trunk)
  • Revision ID: corrado.maurini@upmc.fr-20121218121608-nk82ly9jgsld9u84
updating with trunk, fix uint in TAO solver and hacking the check for tao FindTAO.cmake

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/generation/CSGCGALMeshGenerator2D.cpp
 
1
// Copyright (C) 2012 Johannes Ring
 
2
//
 
3
// This file is part of DOLFIN.
 
4
//
 
5
// DOLFIN is free software: you can redistribute it and/or modify
 
6
// it under the terms of the GNU Lesser General Public License as published by
 
7
// the Free Software Foundation, either version 3 of the License, or
 
8
// (at your option) any later version.
 
9
//
 
10
// DOLFIN is distributed in the hope that it will be useful,
 
11
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 
13
// GNU Lesser General Public License for more details.
 
14
//
 
15
// You should have received a copy of the GNU Lesser General Public License
 
16
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 
17
//
 
18
// Modified by Joachim B Haga 2012
 
19
// Modified by Benjamin Kehlet 2012
 
20
//
 
21
// First added:  2012-05-10
 
22
// Last changed: 2012-11-14
 
23
 
 
24
#include <vector>
 
25
#include <cmath>
 
26
 
 
27
#include <dolfin/common/constants.h>
 
28
#include <dolfin/math/basic.h>
 
29
#include <dolfin/mesh/Mesh.h>
 
30
#include <dolfin/mesh/MeshEditor.h>
 
31
 
 
32
#include "CSGCGALMeshGenerator2D.h"
 
33
#include "CSGGeometry.h"
 
34
#include "CSGOperators.h"
 
35
#include "CSGPrimitives2D.h"
 
36
 
 
37
#ifdef HAS_CGAL
 
38
#include "CGALMeshBuilder.h"
 
39
 
 
40
#include <CGAL/basic.h>
 
41
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
42
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
43
 
 
44
#include <CGAL/Gmpq.h>
 
45
#include <CGAL/Lazy_exact_nt.h>
 
46
#include <CGAL/Simple_cartesian.h>
 
47
#include <CGAL/Bounded_kernel.h>
 
48
#include <CGAL/Nef_polyhedron_2.h>
 
49
#include <CGAL/Polygon_2.h>
 
50
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
 
51
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
 
52
#include <CGAL/Triangulation_face_base_with_info_2.h>
 
53
#include <CGAL/Delaunay_mesher_2.h>
 
54
#include <CGAL/Delaunay_mesh_face_base_2.h>
 
55
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
 
56
 
 
57
typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_Kernel;
 
58
typedef CGAL::Exact_predicates_inexact_constructions_kernel Inexact_Kernel;
 
59
 
 
60
typedef CGAL::Lazy_exact_nt<CGAL::Gmpq> FT;
 
61
typedef CGAL::Simple_cartesian<FT> EKernel;
 
62
typedef CGAL::Bounded_kernel<EKernel> Extended_kernel;
 
63
typedef CGAL::Nef_polyhedron_2<Extended_kernel> Nef_polyhedron_2;
 
64
typedef Nef_polyhedron_2::Point Nef_point_2;
 
65
 
 
66
typedef Nef_polyhedron_2::Explorer Explorer;
 
67
typedef Explorer::Face_const_iterator Face_const_iterator;
 
68
typedef Explorer::Hole_const_iterator Hole_const_iterator;
 
69
typedef Explorer::Halfedge_around_face_const_circulator Halfedge_around_face_const_circulator;
 
70
typedef Explorer::Vertex_const_handle Vertex_const_handle;
 
71
typedef Explorer::Halfedge_const_handle Halfedge_const_handle;
 
72
 
 
73
typedef CGAL::Triangulation_vertex_base_2<Inexact_Kernel>  Vertex_base;
 
74
typedef CGAL::Constrained_triangulation_face_base_2<Inexact_Kernel> Face_base;
 
75
 
 
76
template <class Gt,
 
77
          class Fb >
 
78
class Enriched_face_base_2 : public Fb {
 
79
public:
 
80
  typedef Gt Geom_traits;
 
81
  typedef typename Fb::Vertex_handle Vertex_handle;
 
82
  typedef typename Fb::Face_handle Face_handle;
 
83
 
 
84
  template < typename TDS2 >
 
85
  struct Rebind_TDS {
 
86
    typedef typename Fb::template Rebind_TDS<TDS2>::Other Fb2;
 
87
    typedef Enriched_face_base_2<Gt,Fb2> Other;
 
88
  };
 
89
 
 
90
protected:
 
91
  int status;
 
92
 
 
93
public:
 
94
  Enriched_face_base_2(): Fb(), status(-1) {};
 
95
 
 
96
  Enriched_face_base_2(Vertex_handle v0,
 
97
                       Vertex_handle v1,
 
98
                       Vertex_handle v2)
 
99
    : Fb(v0,v1,v2), status(-1) {};
 
100
 
 
101
  Enriched_face_base_2(Vertex_handle v0,
 
102
                       Vertex_handle v1,
 
103
                       Vertex_handle v2,
 
104
                       Face_handle n0,
 
105
                       Face_handle n1,
 
106
                       Face_handle n2)
 
107
    : Fb(v0,v1,v2,n0,n1,n2), status(-1) {};
 
108
 
 
109
  inline
 
110
  bool is_in_domain() const { return (status%2 == 1); };
 
111
 
 
112
  inline
 
113
  void set_in_domain(const bool b) { status = (b ? 1 : 0); };
 
114
 
 
115
  inline
 
116
  void set_counter(int i) { status = i; };
 
117
 
 
118
  inline
 
119
  int counter() const { return status; };
 
120
 
 
121
  inline
 
122
  int& counter() { return status; };
 
123
}; // end class Enriched_face_base_2
 
124
 
 
125
typedef CGAL::Triangulation_vertex_base_2<Inexact_Kernel> Vb;
 
126
typedef CGAL::Triangulation_vertex_base_with_info_2<std::size_t, Inexact_Kernel, Vb> Vbb;
 
127
typedef Enriched_face_base_2<Inexact_Kernel, Face_base> Fb;
 
128
typedef CGAL::Triangulation_data_structure_2<Vbb, Fb> TDS;
 
129
typedef CGAL::Exact_predicates_tag Itag;
 
130
typedef CGAL::Constrained_Delaunay_triangulation_2<Inexact_Kernel, TDS, Itag> CDT;
 
131
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Mesh_criteria_2;
 
132
typedef CGAL::Delaunay_mesher_2<CDT, Mesh_criteria_2> CGAL_Mesher_2;
 
133
 
 
134
typedef CDT::Vertex_handle Vertex_handle;
 
135
typedef CDT::Face_handle Face_handle;
 
136
typedef CDT::All_faces_iterator All_faces_iterator;
 
137
 
 
138
typedef CGAL::Polygon_2<Inexact_Kernel> Polygon_2;
 
139
typedef Inexact_Kernel::Point_2 Point_2;
 
140
 
 
141
using namespace dolfin;
 
142
 
 
143
//-----------------------------------------------------------------------------
 
144
CSGCGALMeshGenerator2D::CSGCGALMeshGenerator2D(const CSGGeometry& geometry)
 
145
  : geometry(geometry)
 
146
{
 
147
  parameters = default_parameters();
 
148
}
 
149
//-----------------------------------------------------------------------------
 
150
CSGCGALMeshGenerator2D::~CSGCGALMeshGenerator2D() {}
 
151
//-----------------------------------------------------------------------------
 
152
Nef_polyhedron_2 make_circle(const Circle* c)
 
153
{
 
154
  std::vector<Nef_point_2> pts;
 
155
 
 
156
  for (std::size_t i = 0; i < c->fragments(); i++)
 
157
  {
 
158
    const double phi = (2*DOLFIN_PI*i) / c->fragments();
 
159
    const double x = c->center().x() + c->radius()*cos(phi);
 
160
    const double y = c->center().y() + c->radius()*sin(phi);
 
161
    pts.push_back(Nef_point_2(x, y));
 
162
  }
 
163
 
 
164
  return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
 
165
}
 
166
//-----------------------------------------------------------------------------
 
167
Nef_polyhedron_2 make_ellipse(const Ellipse* e)
 
168
{
 
169
  std::vector<Nef_point_2> pts;
 
170
 
 
171
  for (std::size_t i = 0; i < e->fragments(); i++)
 
172
  {
 
173
    const double phi = (2*DOLFIN_PI*i) / e->fragments();
 
174
    const double x = e->center().x() + e->a()*cos(phi);
 
175
    const double y = e->center().y() + e->b()*sin(phi);
 
176
    pts.push_back(Nef_point_2(x, y));
 
177
  }
 
178
 
 
179
  return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
 
180
}
 
181
//-----------------------------------------------------------------------------
 
182
Nef_polyhedron_2 make_rectangle(const Rectangle* r)
 
183
{
 
184
  const double x0 = std::min(r->first_corner().x(), r->first_corner().y());
 
185
  const double y0 = std::max(r->first_corner().x(), r->first_corner().y());
 
186
 
 
187
  const double x1 = std::min(r->second_corner().x(), r->second_corner().y());
 
188
  const double y1 = std::max(r->second_corner().x(), r->second_corner().y());
 
189
 
 
190
  std::vector<Nef_point_2> pts;
 
191
  pts.push_back(Nef_point_2(x0, x1));
 
192
  pts.push_back(Nef_point_2(y0, x1));
 
193
  pts.push_back(Nef_point_2(y0, y1));
 
194
  pts.push_back(Nef_point_2(x0, y1));
 
195
 
 
196
  return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
 
197
}
 
198
//-----------------------------------------------------------------------------
 
199
Nef_polyhedron_2 make_polygon(const Polygon* p)
 
200
{
 
201
  std::vector<Nef_point_2> pts;
 
202
  std::vector<Point>::const_iterator v;
 
203
  for (v = p->vertices().begin(); v != p->vertices().end(); ++v)
 
204
    pts.push_back(Nef_point_2(v->x(), v->y()));
 
205
 
 
206
  return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
 
207
}
 
208
//-----------------------------------------------------------------------------
 
209
static Nef_polyhedron_2 convertSubTree(const CSGGeometry *geometry)
 
210
{
 
211
  switch (geometry->getType()) {
 
212
  case CSGGeometry::Union:
 
213
  {
 
214
    const CSGUnion* u = dynamic_cast<const CSGUnion*>(geometry);
 
215
    dolfin_assert(u);
 
216
    return convertSubTree(u->_g0.get()) + convertSubTree(u->_g1.get());
 
217
    break;
 
218
  }
 
219
  case CSGGeometry::Intersection:
 
220
  {
 
221
    const CSGIntersection* u = dynamic_cast<const CSGIntersection*>(geometry);
 
222
    dolfin_assert(u);
 
223
    return convertSubTree(u->_g0.get()) * convertSubTree(u->_g1.get());
 
224
    break;
 
225
  }
 
226
  case CSGGeometry::Difference:
 
227
  {
 
228
    const CSGDifference* u = dynamic_cast<const CSGDifference*>(geometry);
 
229
    dolfin_assert(u);
 
230
    return convertSubTree(u->_g0.get()) - convertSubTree(u->_g1.get());
 
231
    break;
 
232
  }
 
233
  case CSGGeometry::Circle:
 
234
  {
 
235
    const Circle* c = dynamic_cast<const Circle*>(geometry);
 
236
    dolfin_assert(c);
 
237
    return make_circle(c);
 
238
    break;
 
239
  }
 
240
  case CSGGeometry::Ellipse:
 
241
  {
 
242
    const Ellipse* c = dynamic_cast<const Ellipse*>(geometry);
 
243
    dolfin_assert(c);
 
244
    return make_ellipse(c);
 
245
    break;
 
246
  }
 
247
  case CSGGeometry::Rectangle:
 
248
  {
 
249
    const Rectangle* r = dynamic_cast<const Rectangle*>(geometry);
 
250
    dolfin_assert(r);
 
251
    return make_rectangle(r);
 
252
    break;
 
253
  }
 
254
  case CSGGeometry::Polygon:
 
255
  {
 
256
    const Polygon* p = dynamic_cast<const Polygon*>(geometry);
 
257
    dolfin_assert(p);
 
258
    return make_polygon(p);
 
259
    break;
 
260
  }
 
261
  default:
 
262
    dolfin_error("CSGCGALMeshGenerator2D.cpp",
 
263
                 "converting geometry to cgal polyhedron",
 
264
                 "Unhandled primitive type");
 
265
  }
 
266
  return Nef_polyhedron_2();
 
267
}
 
268
//-----------------------------------------------------------------------------
 
269
// Taken from examples/Triangulation_2/polygon_triangulation.cpp in the
 
270
// CGAL source tree.
 
271
//
 
272
// Explore set of facets connected with non constrained edges,
 
273
// and attribute to each such set a counter.
 
274
// We start from facets incident to the infinite vertex, with a counter
 
275
// set to 0. Then we recursively consider the non-explored facets incident
 
276
// to constrained edges bounding the former set and increase thecounter by 1.
 
277
// Facets in the domain are those with an odd nesting level.
 
278
void mark_domains(CDT& ct,
 
279
                  CDT::Face_handle start,
 
280
                  int index,
 
281
                  std::list<CDT::Edge>& border)
 
282
{
 
283
  if (start->counter() != -1)
 
284
    return;
 
285
 
 
286
  std::list<CDT::Face_handle> queue;
 
287
  queue.push_back(start);
 
288
 
 
289
  while (!queue.empty())
 
290
  {
 
291
    CDT::Face_handle fh = queue.front();
 
292
    queue.pop_front();
 
293
    if (fh->counter() == -1)
 
294
    {
 
295
      fh->counter() = index;
 
296
      fh->set_in_domain(index%2 == 1);
 
297
      for (int i = 0; i < 3; i++)
 
298
      {
 
299
        CDT::Edge e(fh, i);
 
300
        CDT::Face_handle n = fh->neighbor(i);
 
301
        if (n->counter() == -1)
 
302
        {
 
303
          if (ct.is_constrained(e))
 
304
            border.push_back(e);
 
305
          else
 
306
            queue.push_back(n);
 
307
        }
 
308
      }
 
309
    }
 
310
  }
 
311
}
 
312
//-----------------------------------------------------------------------------
 
313
void mark_domains(CDT& cdt)
 
314
{
 
315
  for (CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it)
 
316
    it->set_counter(-1);
 
317
 
 
318
  int index = 0;
 
319
  std::list<CDT::Edge> border;
 
320
  mark_domains(cdt, cdt.infinite_face(), index++, border);
 
321
  while (!border.empty())
 
322
  {
 
323
    CDT::Edge e = border.front();
 
324
    border.pop_front();
 
325
    CDT::Face_handle n = e.first->neighbor(e.second);
 
326
    if (n->counter() == -1)
 
327
      mark_domains(cdt, n, e.first->counter()+1, border);
 
328
  }
 
329
}
 
330
//-----------------------------------------------------------------------------
 
331
void CSGCGALMeshGenerator2D::generate(Mesh& mesh)
 
332
{
 
333
  Nef_polyhedron_2 cgal_geometry = convertSubTree(&geometry);
 
334
 
 
335
  // Create empty CGAL triangulation
 
336
  CDT cdt;
 
337
 
 
338
  // Explore the Nef polyhedron and insert constraints in the triangulation
 
339
  Explorer explorer = cgal_geometry.explorer();
 
340
  Face_const_iterator fit = explorer.faces_begin();
 
341
  for (; fit != explorer.faces_end(); fit++)
 
342
  {
 
343
    // Skip face if it is not part of polygon
 
344
    if (!explorer.mark(fit))
 
345
      continue;
 
346
 
 
347
    Halfedge_around_face_const_circulator hafc = explorer.face_cycle(fit), done(hafc);
 
348
    do {
 
349
      Vertex_handle va = cdt.insert(Point_2(to_double(hafc->vertex()->point().x()),
 
350
                                            to_double(hafc->vertex()->point().y())));
 
351
      Vertex_handle vb = cdt.insert(Point_2(to_double(hafc->next()->vertex()->point().x()),
 
352
                                            to_double(hafc->next()->vertex()->point().y())));
 
353
      cdt.insert_constraint(va, vb);
 
354
      hafc++;
 
355
    } while (hafc != done);
 
356
 
 
357
    Hole_const_iterator hit = explorer.holes_begin(fit);
 
358
    for (; hit != explorer.holes_end(fit); hit++)
 
359
    {
 
360
      Halfedge_around_face_const_circulator hafc(hit), done(hit);
 
361
      do {
 
362
        Vertex_handle va = cdt.insert(Point_2(to_double(hafc->vertex()->point().x()),
 
363
                                              to_double(hafc->vertex()->point().y())));
 
364
        Vertex_handle vb = cdt.insert(Point_2(to_double(hafc->next()->vertex()->point().x()),
 
365
                                              to_double(hafc->next()->vertex()->point().y())));
 
366
        cdt.insert_constraint(va, vb);
 
367
        hafc++;
 
368
      } while (hafc != done);
 
369
    }
 
370
  }
 
371
 
 
372
  // Mark parts that are inside and outside the domain
 
373
  mark_domains(cdt);
 
374
 
 
375
  // Create mesher
 
376
  CGAL_Mesher_2 mesher(cdt);
 
377
 
 
378
  // Add seeds for all faces in the domain
 
379
  std::list<Point_2> list_of_seeds;
 
380
  for(CDT::Finite_faces_iterator fit = cdt.finite_faces_begin();
 
381
      fit != cdt.finite_faces_end(); ++fit)
 
382
  {
 
383
    if (fit->is_in_domain())
 
384
    {
 
385
      // Calculate center of triangle and add to list of seeds
 
386
      Point_2 p0 = fit->vertex(0)->point();
 
387
      Point_2 p1 = fit->vertex(1)->point();
 
388
      Point_2 p2 = fit->vertex(2)->point();
 
389
      const double x = (p0[0] + p1[0] + p2[0]) / 3;
 
390
      const double y = (p0[1] + p1[1] + p2[1]) / 3;
 
391
      list_of_seeds.push_back(Point_2(x, y));
 
392
    }
 
393
  }
 
394
  mesher.set_seeds(list_of_seeds.begin(), list_of_seeds.end(), true);
 
395
 
 
396
  // Set shape and size criteria
 
397
  Mesh_criteria_2 criteria(parameters["triangle_shape_bound"],
 
398
                           parameters["cell_size"]);
 
399
  mesher.set_criteria(criteria);
 
400
 
 
401
  // Refine CGAL mesh/triangulation
 
402
  mesher.refine_mesh();
 
403
 
 
404
  // Make sure triangulation is valid
 
405
  dolfin_assert(cdt.is_valid());
 
406
 
 
407
  // Clear mesh
 
408
  mesh.clear();
 
409
 
 
410
  // Get various dimensions
 
411
  const std::size_t gdim = cdt.finite_vertices_begin()->point().dimension();
 
412
  const std::size_t tdim = cdt.dimension();
 
413
  const std::size_t num_vertices = cdt.number_of_vertices();
 
414
 
 
415
  // Count valid cells
 
416
  std::size_t num_cells = 0;
 
417
  CDT::Finite_faces_iterator cgal_cell;
 
418
  for (cgal_cell = cdt.finite_faces_begin(); cgal_cell != cdt.finite_faces_end(); ++cgal_cell)
 
419
  {
 
420
    // Add cell if it is in the domain
 
421
    if (cgal_cell->is_in_domain())
 
422
    {
 
423
      num_cells++;
 
424
    }
 
425
  }
 
426
 
 
427
  // Create a MeshEditor and open
 
428
  dolfin::MeshEditor mesh_editor;
 
429
  mesh_editor.open(mesh, tdim, gdim);
 
430
  mesh_editor.init_vertices(num_vertices);
 
431
  mesh_editor.init_cells(num_cells);
 
432
 
 
433
  // Add vertices to mesh
 
434
  std::size_t vertex_index = 0;
 
435
  CDT::Finite_vertices_iterator cgal_vertex;
 
436
  for (cgal_vertex = cdt.finite_vertices_begin();
 
437
       cgal_vertex != cdt.finite_vertices_end(); ++cgal_vertex)
 
438
  {
 
439
    // Get vertex coordinates and add vertex to the mesh
 
440
    Point p;
 
441
    p[0] = cgal_vertex->point()[0];
 
442
    p[1] = cgal_vertex->point()[1];
 
443
    if (gdim == 3)
 
444
      p[2] = cgal_vertex->point()[2];
 
445
 
 
446
    // Add mesh vertex
 
447
    mesh_editor.add_vertex(vertex_index, p);
 
448
 
 
449
    // Attach index to vertex and increment
 
450
    cgal_vertex->info() = vertex_index++;
 
451
  }
 
452
  dolfin_assert(vertex_index == num_vertices);
 
453
 
 
454
  // Add cells to mesh
 
455
  std::size_t cell_index = 0;
 
456
  for (cgal_cell = cdt.finite_faces_begin(); cgal_cell != cdt.finite_faces_end(); ++cgal_cell)
 
457
  {
 
458
    // Add cell if it is in the domain
 
459
    if (cgal_cell->is_in_domain())
 
460
    {
 
461
      mesh_editor.add_cell(cell_index++,
 
462
                           cgal_cell->vertex(0)->info(),
 
463
                           cgal_cell->vertex(1)->info(),
 
464
                           cgal_cell->vertex(2)->info());
 
465
    }
 
466
  }
 
467
  dolfin_assert(cell_index == num_cells);
 
468
 
 
469
  // Close mesh editor
 
470
  mesh_editor.close();
 
471
 
 
472
  // Build DOLFIN mesh from CGAL triangulation
 
473
  // FIXME: Why does this not work correctly?
 
474
  //CGALMeshBuilder::build(mesh, cdt);
 
475
}
 
476
 
 
477
#else
 
478
namespace dolfin
 
479
{
 
480
  CSGCGALMeshGenerator2D::CSGCGALMeshGenerator2D(const CSGGeometry& geometry)
 
481
    :geometry(geometry)
 
482
  {
 
483
    dolfin_error("CSGCGALMeshGenerator2D.cpp",
 
484
                 "Create mesh generator",
 
485
                 "Dolfin must be compiled with CGAL to use this feature.");
 
486
  }
 
487
  //-----------------------------------------------------------------------------
 
488
  CSGCGALMeshGenerator2D::~CSGCGALMeshGenerator2D(){}
 
489
  //-----------------------------------------------------------------------------
 
490
  void CSGCGALMeshGenerator2D::generate(Mesh& mesh) {}
 
491
}
 
492
 
 
493
 
 
494
#endif
 
495
//-----------------------------------------------------------------------------