~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to demo/mesh/intersection/cpp/main.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
ImportĀ upstreamĀ versionĀ 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
// Copyright (C) 2008 Kristoffer Selim.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2008-10-08
 
5
// Last changed: 2008-12-18
 
6
//
 
7
// Modified by Anders Logg, 2008.
 
8
 
 
9
#include <dolfin.h>
 
10
#include <math.h>
 
11
 
 
12
using namespace dolfin;
 
13
 
 
14
int main()
 
15
{
 
16
  // Create meshes (omega0 overlapped by omega1)
 
17
  UnitCircle omega0(20);
 
18
  UnitSquare omega1(20, 20);
 
19
    
 
20
  // Access mesh geometry
 
21
  MeshGeometry& geometry = omega0.geometry();
 
22
 
 
23
  // Move and scale circle
 
24
  for (VertexIterator vertex(omega0); !vertex.end(); ++vertex)
 
25
  {
 
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;
 
29
  }
 
30
 
 
31
  // Iterate over angle
 
32
  const double dtheta = 0.1*DOLFIN_PI;
 
33
  for (double theta = 0; theta < 2*DOLFIN_PI; theta += dtheta)
 
34
  {
 
35
    // Compute intersection with boundary of square
 
36
    BoundaryMesh boundary(omega1);
 
37
    std::vector<unsigned int> cells;
 
38
    omega0.intersection(boundary, cells, false);
 
39
    
 
40
    // Copy values to mesh function for plotting
 
41
    MeshFunction<unsigned int> intersection(omega0, omega0.topology().dim());
 
42
    intersection = 0;
 
43
    for (unsigned int i = 0; i < cells.size(); i++)
 
44
      intersection.set(cells[i], 1);
 
45
    
 
46
    // Plot intersection
 
47
    plot(intersection);
 
48
 
 
49
    // Rotate circle around (0.5, 0.5)
 
50
    for (VertexIterator vertex(omega0); !vertex.end(); ++vertex)
 
51
    {
 
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);
 
57
    }
 
58
  }
 
59
}