1
// Copyright (C) 2008 Kristoffer Selim.
2
// Licensed under the GNU LGPL Version 2.1.
4
// First added: 2008-10-08
5
// Last changed: 2008-12-18
7
// Modified by Anders Logg, 2008.
12
using namespace dolfin;
16
// Create meshes (omega0 overlapped by omega1)
17
UnitCircle omega0(20);
18
UnitSquare omega1(20, 20);
20
// Access mesh geometry
21
MeshGeometry& geometry = omega0.geometry();
23
// Move and scale circle
24
for (VertexIterator vertex(omega0); !vertex.end(); ++vertex)
26
double* x = geometry.x(vertex->index());
27
x[0] = 0.5*x[0] + 1.0;
28
x[1] = 0.5*x[1] + 1.0;
32
const double dtheta = 0.1*DOLFIN_PI;
33
for (double theta = 0; theta < 2*DOLFIN_PI; theta += dtheta)
35
// Compute intersection with boundary of square
36
BoundaryMesh boundary(omega1);
37
std::vector<unsigned int> cells;
38
omega0.intersection(boundary, cells, false);
40
// Copy values to mesh function for plotting
41
MeshFunction<unsigned int> intersection(omega0, omega0.topology().dim());
43
for (unsigned int i = 0; i < cells.size(); i++)
44
intersection.set(cells[i], 1);
49
// Rotate circle around (0.5, 0.5)
50
for (VertexIterator vertex(omega0); !vertex.end(); ++vertex)
52
double* x = geometry.x(vertex->index());
53
const double xr = x[0] - 0.5;
54
const double yr = x[1] - 0.5;
55
x[0] = 0.5 + (cos(dtheta)*xr - sin(dtheta)*yr);
56
x[1] = 0.5 + (sin(dtheta)*xr + cos(dtheta)*yr);