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; }Handbook
Last modified: Wed Jun 6 10:14:54 PDT 2001