← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to dolfin/la/SingularSolver.h

  • Committer: Anders Logg
  • Date: 2008-05-15 16:02:44 UTC
  • Revision ID: logg@simula.no-20080515160244-e2a4wkfxfzmts16s
Add new class SingularSolver:

This class provides a linear solver for singular linear systems
Ax = b where A has a one-dimensional null-space (kernel). This
may happen for example when solving Poisson's equation with
pure Neumann boundary conditions.

The solver attempts to create an extended non-singular system
by adding the constraint [1, 1, 1, ...]^T x = 0.

If an optional mass matrix M is supplied, the solver attempts
to create an extended non-singular system by adding the
constraint m^T x = 0 where m is the lumped mass matrix. This
corresponds to setting the average (integral) of the finite
element function with coefficients x to zero.

The solver makes not attempt to check that the null-space is
indeed one-dimensional. It is also assumed that the system
Ax = b retains its sparsity pattern between calls to solve().

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/la/SingularSolver.h
 
1
// Copyright (C) 2008 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2005-09-19
 
5
// Last changed: 2008-05-10
 
6
 
 
7
#ifndef __SINGULAR_SOLVER_H
 
8
#define __SINGULAR_SOLVER_H
 
9
 
 
10
#include <dolfin/common/types.h>
 
11
#include <dolfin/parameter/Parametrized.h>
 
12
#include "LinearSolver.h"
 
13
 
 
14
namespace dolfin
 
15
{
 
16
  /// This class provides a linear solver for singular linear systems
 
17
  /// Ax = b where A has a one-dimensional null-space (kernel). This
 
18
  /// may happen for example when solving Poisson's equation with
 
19
  /// pure Neumann boundary conditions.
 
20
  ///
 
21
  /// The solver attempts to create an extended non-singular system
 
22
  /// by adding the constraint [1, 1, 1, ...]^T x = 0.
 
23
  ///
 
24
  /// If an optional mass matrix M is supplied, the solver attempts
 
25
  /// to create an extended non-singular system by adding the
 
26
  /// constraint m^T x = 0 where m is the lumped mass matrix. This
 
27
  /// corresponds to setting the average (integral) of the finite
 
28
  /// element function with coefficients x to zero.
 
29
  ///
 
30
  /// The solver makes not attempt to check that the null-space is
 
31
  /// indeed one-dimensional. It is also assumed that the system
 
32
  /// Ax = b retains its sparsity pattern between calls to solve().
 
33
  
 
34
  class SingularSolver : public Parametrized
 
35
  {
 
36
  public:
 
37
 
 
38
    /// Create linear solver
 
39
    SingularSolver(SolverType solver_type=lu, PreconditionerType pc_type=ilu);
 
40
 
 
41
    /// Destructor
 
42
    ~SingularSolver();
 
43
 
 
44
    /// Solve linear system Ax = b
 
45
    uint solve(const GenericMatrix& A, GenericVector& x, const GenericVector& b);
 
46
 
 
47
    /// Solve linear system Ax = b using mass matrix M for setting constraint
 
48
    uint solve(const GenericMatrix& A, GenericVector& x, const GenericVector& b,
 
49
               const GenericMatrix& M);
 
50
 
 
51
  private:
 
52
 
 
53
    // Initialize extended system
 
54
    void init(const GenericMatrix& A);
 
55
 
 
56
    // Create extended system
 
57
    void create(const GenericMatrix& A, const GenericVector& b, const GenericMatrix* M);
 
58
 
 
59
    // Linear solver
 
60
    LinearSolver linear_solver;
 
61
 
 
62
    // Extended matrix
 
63
    GenericMatrix* B;
 
64
 
 
65
    // Solution of extended system
 
66
    GenericVector* y;
 
67
    
 
68
    // Extended vector
 
69
    GenericVector* c;
 
70
 
 
71
  };
 
72
 
 
73
}
 
74
 
 
75
#endif