Simple Example

The code below shows a simple code for solving Poisson's equation by cg-iteration. The pure code can be found here .

#include "source/ugblock.h"

using namespace ::_COLSAMM_;



          /* Poisson on a Ball */


int main()
{
  //  cout.precision(10);
  //  cout.setf(ios::fixed,ios::floatfield);
  int iteration, i, n;
  double delta, delta_prime, beta, tau, eps, h;              

  n = 10;                                      // number of grid points
  iteration = 100;
  eps = 0.0000000001;

  cout << "\n Solving Poisson's equation \n";
  cout << " with Dirichlet boundary conditions on a cylinder by cg: ";
  cout << "\n Solving Poisson's equation on a cylinder by cg: ";
  cout << "\n Meshsize of grid: " << h;
  cout << "\n Maximal number of iterations: " << iteration;
  cout << "\n Stopping parameter for cg: " << eps << endl;

  // Part 1: Construction of a domain and markers:
  //-----------------------------------------------
  Cylinder cylinder(2.0,1.0,4.0);

  Boundary_Marker boundary(&cylinder);
  Unstructured_grid_Marker inner_points(&cylinder);
  inner_points.complement(&boundary);


  // Part 2: Construction of a grid:
  //--------------------------------
  Blockgrid grid(&cylinder,n);                    // Construction of grid

  // Part 3: Definition Variables and other object which need storage:
  //------------------------------------------------------------------
  Variable<double> f(grid);
  Variable<double> g(grid);
  Variable<double> d(grid);               

  Variable<double> r(grid);
  Variable<double> e(grid);
  Variable<double> u(grid);
  Variable<double> u_exakt(grid);


  // Part 4: The mathematical code:
  //-------------------------------
  Function1<double> Sin(sin);                   // definition of functions
  Function1<double> Sinh(sinh);

  X_coordinate X(grid);
  Y_coordinate Y(grid);
  Z_coordinate Z(grid);

  Local_stiffness_matrix<double> Laplace(grid);
  Local_stiffness_matrix<double> Helm(grid);
                  
  Laplace.Calculate(grad(v_())*grad(w_()));
  Helm.Calculate(v_()*w_());

  u_exakt = Sin(X)*Sinh(Y)*Z;     // exact solution
  
  u = 0.0;
  u = u_exakt | boundary;         // setting Dirichlet boundary values
  
  r = 0.0;                        // right hand side is zero
  f = Helm(f)                          | inner_points;
  
  cout << "\n cg: " << endl;
  cout << "-------------------- " << endl;
  
  r = Laplace(u) - f                   | inner_points;
  d = -1.0 * r                         | inner_points;
  
  delta = product(r,r,inner_points);
  
  for(i=1;i<=iteration && delta > eps;++i) {
    cout << i << ". ";
    
    g = Laplace(d)                     | inner_points;
    
    tau = delta / product(d,g,inner_points);
    
    r = r + tau * g                    | inner_points;
    u = u + tau * d                    | inner_points;

    delta_prime = product(r,r,inner_points);
    beta = delta_prime / delta;   
    delta = delta_prime;
    d = beta*d - r                     | inner_points;
    e = u - u_exakt;

    cout << "Iteration: Error: Max:" << L_infty(e) << endl;
  }

  cout << "\n end. " << endl;
}                                                 

Last modified: Wed Jun 6 10:14:54 PDT 2001