1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
// Jacobian.hh -- data structures for the Jacobian matrix
// $Id$
//
// Jacobian -- abstract base class to describe a Jacobian matrix
//
// dense_Jacobian - Jacobian stored as a dense matrix
//
//
// prerequisites:
// "stdc.h"
// "config.hh"
// "../jtutil/util.hh"
// "../jtutil/array.hh"
// "../jtutil/cpm_map.hh"
// "../jtutil/linear_map.hh"
// "coords.hh"
// "grid.hh"
// "fd_grid.hh"
// "patch.hh"
// "patch_edge.hh"
// "patch_interp.hh"
// "ghost_zone.hh"
// "patch_system.hh"
//
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// A Jacobian object stores the (a) Jacobian matrix for a patch system.
// Jacobian is an abstract base class, from which we derive a separate
// concrete class for each type of sparsity pattern (more precisely for
// each type of storage scheme) of the Jacobian matrix.
//
// A row/column index of the Jacobian (denoted II/JJ) is a 0-origin grid
// point number within the patch system.
//
// Note that the APIs here implicitly assume there is only a *single* gridfn
// in the Jacobian computation. (If we had N gridfns for this, then the
// Jacobian would really be a block matrix with N*N blocks.) This isn't
// very general, but matches our present use in this thorn.
//
class Jacobian
{
public:
// basic meta-info
patch_system& my_patch_system();
public:
// access a matrix element via row/column indices
virtual const fp& operator()(int II, int JJ) const = 0; // rvalue
virtual fp& operator()(int II, int JJ) = 0; // lvalue
// access a matrix element via row index and column (patch,irho,isigma)
const fp& operator() // rvalue
(int II, const patch& JJ_patch, int JJ_irho, int JJ_isigma)
const
{
const int JJ
= ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
return operator()(II,JJ);
}
fp& operator() // lvalue
(int II, const patch& JJ_patch, int JJ_irho, int JJ_isigma)
{
const int JJ
= ps_.gpn_of_patch_irho_isigma(JJ_patch, JJ_irho,JJ_isigma);
return operator()(II,JJ);
}
// zero the whole matrix or a single row
virtual void zero() = 0;
virtual void zero_row(int II) = 0;
// basic meta-info
const patch_system& my_patch_system() const { return ps_; }
// constructor, destructor
Jacobian(const patch_system& ps);
virtual ~Jacobian() { }
private:
const patch_system& ps_;
};
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This class stores the Jacobian as a dense matrix in Fortran (column)
// order.
//
class dense_Jacobian
: public Jacobian
{
public:
// access a matrix element via row/column indices
const fp& operator()(int II, int JJ) const // rvalue
{ return matrix_(JJ,II); }
fp& operator()(int II, int JJ) // lvalue
{ return matrix_(JJ,II); }
// zero the whole matrix or a single row
void zero();
void zero_row(int II);
// constructor, destructor
public:
dense_Jacobian(const patch_system& ps);
~dense_Jacobian() { }
private:
// subscripts are (JJ,II)
jtutil::array2d<fp> matrix_;
};
|
