← Back to branch summary

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

« back to all changes in this revision

Viewing changes to demo/pde/stokes/taylor-hood/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:
demo/pde/stokes/taylor-hood/cpp/main.cpp
 
1
// Copyright (C) 2006-2008 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2006-02-09
 
5
// Last changed: 2008-12-12
 
6
//
 
7
// This demo solves the Stokes equations, using quadratic elements for
 
8
// the velocity and first degree elements for the pressure
 
9
// (Taylor-Hood elements). The sub domains for the different boundary
 
10
// conditions used in this simulation are computed by the demo program
 
11
// in src/demo/mesh/subdomains.
 
12
 
 
13
#include <dolfin.h>
 
14
#include "Stokes.h"
 
15
 
 
16
using namespace dolfin;
 
17
 
 
18
int main()
 
19
{
 
20
  // Function for no-slip boundary condition for velocity
 
21
  class Noslip : public Function
 
22
  {
 
23
    void eval(double* values, const double* x) const
 
24
    {
 
25
      values[0] = 0.0;
 
26
      values[1] = 0.0;
 
27
    }
 
28
  };
 
29
 
 
30
  // Function for inflow boundary condition for velocity
 
31
  class Inflow : public Function
 
32
  {
 
33
    void eval(double* values, const double* x) const
 
34
    {
 
35
      values[0] = -sin(x[1]*DOLFIN_PI);
 
36
      values[1] = 0.0;
 
37
    }
 
38
  };
 
39
 
 
40
  // Read mesh and sub domain markers
 
41
  Mesh mesh("../../../../../data/meshes/dolfin-2.xml.gz");
 
42
  mesh.init(1);
 
43
  mesh.order();
 
44
  MeshFunction<unsigned int> sub_domains(mesh, "../subdomains.xml.gz");
 
45
 
 
46
  // Create function space and subspaces
 
47
  StokesFunctionSpace W(mesh);
 
48
  SubSpace W0(W, 0);
 
49
  SubSpace W1(W, 1);
 
50
 
 
51
  // Create functions for boundary conditions
 
52
  Noslip noslip;
 
53
  Inflow inflow;
 
54
  Constant zero(0.0);
 
55
  
 
56
  // No-slip boundary condition for velocity
 
57
  DirichletBC bc0(W0, noslip, sub_domains, 0);
 
58
 
 
59
  // Inflow boundary condition for velocity
 
60
  DirichletBC bc1(W0, inflow, sub_domains, 1);
 
61
 
 
62
  // Boundary condition for pressure at outflow
 
63
  DirichletBC bc2(W1, zero, sub_domains, 2);
 
64
 
 
65
  // Collect boundary conditions
 
66
  std::vector<BoundaryCondition*> bcs;
 
67
  bcs.push_back(&bc0); bcs.push_back(&bc1); bcs.push_back(&bc2);
 
68
 
 
69
  // Set up PDE
 
70
  Constant f(2, 0.0);
 
71
  StokesBilinearForm a(W, W);
 
72
  StokesLinearForm L(W);
 
73
  L.f = f;
 
74
  VariationalProblem pde(a, L, bcs);
 
75
 
 
76
  // Solve PDE
 
77
  Function w;
 
78
  pde.set("linear solver", "direct");
 
79
  pde.solve(w);
 
80
  Function u = w[0];
 
81
  Function p = w[1];
 
82
 
 
83
  // Plot solution
 
84
  plot(u);
 
85
  plot(p);
 
86
 
 
87
  // Save solution in VTK format
 
88
  File ufile_pvd("velocity.pvd");
 
89
  ufile_pvd << u;
 
90
  File pfile_pvd("pressure.pvd");
 
91
  pfile_pvd << p;
 
92
}