← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to demo/pde/periodic/cpp/main.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/pde/periodic/cpp/main.cpp
 
1
// Copyright (C) 2007 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2007-07-11
 
5
// Last changed: 2007-08-20
 
6
//
 
7
// This demo program solves Poisson's equation
 
8
//
 
9
//     - div grad u(x, y) = f(x, y)
 
10
//
 
11
// on the unit square with homogeneous Dirichlet boundary conditions
 
12
// at y = 0, 1 and periodic boundary conditions at x = 0, 1.
 
13
 
 
14
#include <dolfin.h>
 
15
#include "Poisson.h"
 
16
  
 
17
using namespace dolfin;
 
18
 
 
19
int main()
 
20
{
 
21
  // Source term
 
22
  class Source : public Function
 
23
  {
 
24
  public:
 
25
    
 
26
    Source(Mesh& mesh) : Function(mesh) {}
 
27
 
 
28
    real eval(const real* x) const
 
29
    {
 
30
      real dx = x[0] - 0.5;
 
31
      real dy = x[1] - 0.5;
 
32
      return x[0]*sin(5.0*DOLFIN_PI*x[1]) + 1.0*exp(-(dx*dx + dy*dy)/0.02);
 
33
    }
 
34
 
 
35
  };
 
36
 
 
37
  // Sub domain for Dirichlet boundary condition
 
38
  class DirichletBoundary : public SubDomain
 
39
  {
 
40
    bool inside(const real* x, bool on_boundary) const
 
41
    {
 
42
      return x[1] < DOLFIN_EPS || x[1] > (1.0 - DOLFIN_EPS) && on_boundary;
 
43
    }
 
44
  };
 
45
 
 
46
  // Sub domain for Periodic boundary condition
 
47
  class PeriodicBoundary : public SubDomain
 
48
  {
 
49
    bool inside(const real* x, bool on_boundary) const
 
50
    {
 
51
      return x[0] < DOLFIN_EPS && x[0] > -DOLFIN_EPS && on_boundary;
 
52
    }
 
53
 
 
54
    void map(const real* x, real* y) const
 
55
    {
 
56
      y[0] = x[0] - 1.0;
 
57
      y[1] = x[1];
 
58
    }
 
59
  };
 
60
 
 
61
  // Create mesh
 
62
  UnitSquare mesh(32, 32);
 
63
 
 
64
  // Create functions
 
65
  Source f(mesh);
 
66
 
 
67
  // Create Dirichlet boundary condition
 
68
  Function u0(mesh, 0.0);
 
69
  DirichletBoundary dirichlet_boundary;
 
70
  DirichletBC bc0(u0, mesh, dirichlet_boundary);
 
71
  
 
72
  // Create periodic boundary condition
 
73
  PeriodicBoundary periodic_boundary;
 
74
  PeriodicBC bc1(mesh, periodic_boundary);
 
75
 
 
76
  // Collect boundary conditions
 
77
  Array<BoundaryCondition*> bcs(&bc0, &bc1);
 
78
 
 
79
  // Define PDE
 
80
  PoissonBilinearForm a;
 
81
  PoissonLinearForm L(f);
 
82
  LinearPDE pde(a, L, mesh, bcs);
 
83
 
 
84
  // Solve PDE
 
85
  Function u;
 
86
  pde.solve(u);
 
87
 
 
88
  // Plot solution
 
89
  plot(u);
 
90
 
 
91
  // Save solution to file
 
92
  File file("poisson.pvd");
 
93
  file << u;
 
94
 
 
95
  return 0;
 
96
}