← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty-proposed

« back to all changes in this revision

Viewing changes to nfem/lib/form_manip.cc

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/lib/form_manip.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
///
 
20
/// =========================================================================
 
21
#include "rheolef/form_manip.h"
 
22
using namespace std;
 
23
 
 
24
form_manip::form_manip()
 
25
  : _verbose(true),
 
26
    _nform(0),
 
27
    _nrowbloc(0),
 
28
    _ncolbloc(0),
 
29
    _a(0)
 
30
{
 
31
}
 
32
form_manip::size_type
 
33
form_manip::nform() const
 
34
{
 
35
    return _nform;
 
36
}
 
37
form_manip::size_type 
 
38
form_manip::nrowbloc() const
 
39
{
 
40
    return _nrowbloc;
 
41
}
 
42
form_manip::size_type 
 
43
form_manip::ncolbloc() const
 
44
{
 
45
    return _ncolbloc;
 
46
}
 
47
form&
 
48
form_manip::bloc(size_type i, size_type j)
 
49
{
 
50
    return _a[_ncolbloc*i+j];
 
51
}
 
52
const form&
 
53
form_manip::bloc (size_type i, size_type j) const
 
54
{
 
55
    return _a[_ncolbloc*i+j];
 
56
}
 
57
form_manip&
 
58
form_manip::operator << (const form& a) 
 
59
{
 
60
    _a.insert (_a.end(),a);
 
61
    _nform++; 
 
62
    return *this;
 
63
}
 
64
form_manip&
 
65
form_manip::operator >> (form& a)
 
66
{
 
67
  if (_nform != _nrowbloc*_ncolbloc) {
 
68
      error_macro("bad bloc number: " << _nform << ", expect " << _nrowbloc*_ncolbloc);
 
69
  }
 
70
  check_macro (_ncolbloc >0 && _nrowbloc > 0, "size(n,m) may be non-zero in form_manip");
 
71
  //
 
72
  // initialize (X,Y) spaces as product spaces
 
73
  //
 
74
  a.X_ = bloc(0,0).X_;
 
75
  for (size_t j = 1; j < _ncolbloc; j++) {
 
76
    a.X_ = a.X_ * bloc(0,j).X_;
 
77
  }
 
78
  a.Y_ = bloc(0,0).Y_;
 
79
  for (size_t i = 1; i < _nrowbloc; i++) {
 
80
    a.Y_ = a.Y_ * bloc(i,0).Y_;
 
81
  }
 
82
  //
 
83
  // big switch: there is not yet a general case !!!
 
84
  //
 
85
  switch(_nrowbloc) {
 
86
 
 
87
  case 1:
 
88
    switch(_ncolbloc) {
 
89
    case 1: // scalaire    -> scalaire
 
90
      a.uu = bloc(0,0).uu ; a.ub = bloc(0,0).ub ;
 
91
      a.bu = bloc(0,0).bu ; a.bb = bloc(0,0).bb ;
 
92
      break ;
 
93
    case 2: // vecteur(2D) -> scalaire
 
94
      a.uu = hcat(bloc(0,0).uu,bloc(0,1).uu) ; a.ub = hcat(bloc(0,0).ub,bloc(0,1).ub) ;
 
95
      a.bu = hcat(bloc(0,0).bu,bloc(0,1).bu) ; a.bb = hcat(bloc(0,0).bb,bloc(0,1).bb) ;
 
96
      break ;
 
97
    case 3: // vecteur(3D) -> scalaire ou tenseur sym(2D) -> scalaire
 
98
      trace_macro("3D not yet: 2d sym. tenseur sym -> scalar");
 
99
      a.uu = hcat(bloc(0,0).uu,hcat(bloc(0,1).uu,bloc(0,2).uu)) ; 
 
100
      a.ub = hcat(bloc(0,0).ub,hcat(bloc(0,1).ub,bloc(0,2).ub)) ;
 
101
      a.bu = hcat(bloc(0,0).bu,hcat(bloc(0,1).bu,bloc(0,2).bu)) ; 
 
102
      a.bb = hcat(bloc(0,0).bb,hcat(bloc(0,1).bb,bloc(0,2).bb)) ;
 
103
      break ;
 
104
    default:
 
105
      error_macro("vector(" << _ncolbloc << "D) is scalar");
 
106
    }
 
107
    break ;
 
108
    
 
109
  case 2:
 
110
    switch(_ncolbloc) {
 
111
    case 1: // scalaire    -> vecteur(2D)
 
112
      a.uu = vcat(bloc(0,0).uu,bloc(1,0).uu) ; a.ub = vcat(bloc(0,0).ub,bloc(1,0).ub) ; 
 
113
      a.bu = vcat(bloc(0,0).bu,bloc(1,0).bu) ; a.bb = vcat(bloc(0,0).bb,bloc(1,0).bb) ; 
 
114
      break ;
 
115
    case 2: // vecteur(2D) -> vecteur(2D) 
 
116
      a.uu = vcat( hcat(bloc(0,0).uu,bloc(0,1).uu) , hcat(bloc(1,0).uu,bloc(1,1).uu) ) ;
 
117
      a.ub = vcat( hcat(bloc(0,0).ub,bloc(0,1).ub) , hcat(bloc(1,0).ub,bloc(1,1).ub) ) ;
 
118
      a.bu = vcat( hcat(bloc(0,0).bu,bloc(0,1).bu) , hcat(bloc(1,0).bu,bloc(1,1).bu) ) ;
 
119
      a.bb = vcat( hcat(bloc(0,0).bb,bloc(0,1).bb) , hcat(bloc(1,0).bb,bloc(1,1).bb) ) ;
 
120
      break ;
 
121
    case 3: // tenseur(2D) -> vecteur(2D)
 
122
      a.uu = vcat( hcat( hcat(bloc(0,0).uu,bloc(0,1).uu) , bloc(0,2).uu ) ,
 
123
                   hcat( hcat(bloc(1,0).uu,bloc(1,1).uu) , bloc(1,2).uu ) );  
 
124
      a.ub = vcat( hcat( hcat(bloc(0,0).ub,bloc(0,1).ub) , bloc(0,2).ub ) ,
 
125
                   hcat( hcat(bloc(1,0).ub,bloc(1,1).ub) , bloc(1,2).ub ) ) ;  
 
126
      a.bu = vcat( hcat( hcat(bloc(0,0).bu,bloc(0,1).bu) , bloc(0,2).bu ) ,
 
127
                   hcat( hcat(bloc(1,0).bu,bloc(1,1).bu) , bloc(1,2).bu ) ) ;  
 
128
      a.bb = vcat( hcat( hcat(bloc(0,0).bb,bloc(0,1).bb) , bloc(0,2).bb ) ,
 
129
                   hcat( hcat(bloc(1,0).bb,bloc(1,1).bb) , bloc(1,2).bb ) ) ;   
 
130
      break ;
 
131
    case 4: // tenseur(2Daxi) -> vecteur(2D)
 
132
      a.uu = vcat( hcat( hcat(bloc(0,0).uu,bloc(0,1).uu) , hcat(bloc(0,2).uu,bloc(0,3).uu) ) ,
 
133
                   hcat( hcat(bloc(1,0).uu,bloc(1,1).uu) , hcat(bloc(1,2).uu,bloc(1,3).uu) ) );  
 
134
      a.ub = vcat( hcat( hcat(bloc(0,0).ub,bloc(0,1).ub) , hcat(bloc(0,2).ub,bloc(0,3).ub) ) ,
 
135
                   hcat( hcat(bloc(1,0).ub,bloc(1,1).ub) , hcat(bloc(1,2).ub,bloc(1,3).ub) ) ) ;  
 
136
      a.bu = vcat( hcat( hcat(bloc(0,0).bu,bloc(0,1).bu) , hcat(bloc(0,2).bu,bloc(0,3).bu) ) ,
 
137
                   hcat( hcat(bloc(1,0).bu,bloc(1,1).bu) , hcat(bloc(1,2).bu,bloc(1,3).bu) ) ) ;  
 
138
      a.bb = vcat( hcat( hcat(bloc(0,0).bb,bloc(0,1).bb) , hcat(bloc(0,2).bb,bloc(0,3).bb) ) ,
 
139
                   hcat( hcat(bloc(1,0).bb,bloc(1,1).bb) , hcat(bloc(1,2).bb,bloc(1,3).bb) ) ) ;   
 
140
      break ;
 
141
    default:
 
142
      error_macro("vector (" << _ncolbloc << "D) is vecteur(2D)");
 
143
    }
 
144
    break ;
 
145
    
 
146
  case 3:
 
147
    switch(_ncolbloc) {
 
148
    case 1: // scalaire -> vecteur(3D) ou scalaire -> tenseur sym(2D)
 
149
      trace_macro("3D not yet: scalaire -> tenseur sym(2D)");
 
150
      a.uu = vcat(bloc(0,0).uu,vcat(bloc(1,0).uu,bloc(2,0).uu)) ; 
 
151
      a.ub = vcat(bloc(0,0).ub,vcat(bloc(1,0).ub,bloc(2,0).ub)) ;
 
152
      a.bu = vcat(bloc(0,0).bu,vcat(bloc(1,0).bu,bloc(2,0).bu)) ; 
 
153
      a.bb = vcat(bloc(0,0).bb,vcat(bloc(1,0).bb,bloc(2,0).bb)) ;      
 
154
      break ;
 
155
    case 2: // vecteur(2D) -> tenseur(2D)
 
156
      a.uu = hcat( vcat( vcat(bloc(0,0).uu,bloc(1,0).uu) , bloc(2,0).uu ) ,
 
157
                   vcat( vcat(bloc(0,1).uu,bloc(1,1).uu) , bloc(2,1).uu ) ) ;  
 
158
      a.ub = hcat( vcat( vcat(bloc(0,0).ub,bloc(1,0).ub) , bloc(2,0).ub ) ,
 
159
                   vcat( vcat(bloc(0,1).ub,bloc(1,1).ub) , bloc(2,1).ub ) ) ;  
 
160
      a.bu = hcat( vcat( vcat(bloc(0,0).bu,bloc(1,0).bu) , bloc(2,0).bu ) ,
 
161
                   vcat( vcat(bloc(0,1).bu,bloc(1,1).bu) , bloc(2,1).bu ) ) ;  
 
162
      a.bb = hcat( vcat( vcat(bloc(0,0).bb,bloc(1,0).bb) , bloc(2,0).bb ) ,
 
163
                   vcat( vcat(bloc(0,1).bb,bloc(1,1).bb) , bloc(2,1).bb ) ) ;  
 
164
      break ;
 
165
    case 3: // tenseur(2D) -> tenseur(2D) et vecteur(3D) -> vecteur(3D)
 
166
      trace_macro("3D in test...");
 
167
      a.uu = hcat( hcat(bloc(0,0).uu,bloc(0,1).uu) , bloc(0,2).uu ) ;
 
168
      a.ub = hcat( hcat(bloc(0,0).ub,bloc(0,1).ub) , bloc(0,2).ub ) ;
 
169
      a.bu = hcat( hcat(bloc(0,0).bu,bloc(0,1).bu) , bloc(0,2).bu ) ;
 
170
      a.bb = hcat( hcat(bloc(0,0).bb,bloc(0,1).bb) , bloc(0,2).bb ) ;
 
171
      a.uu = vcat( a.uu , hcat( hcat(bloc(1,0).uu,bloc(1,1).uu) , bloc(1,2).uu ) ) ;  
 
172
      a.ub = vcat( a.ub , hcat( hcat(bloc(1,0).ub,bloc(1,1).ub) , bloc(1,2).ub ) ) ;  
 
173
      a.bu = vcat( a.bu , hcat( hcat(bloc(1,0).bu,bloc(1,1).bu) , bloc(1,2).bu ) ) ;  
 
174
      a.bb = vcat( a.bb , hcat( hcat(bloc(1,0).bb,bloc(1,1).bb) , bloc(1,2).bb ) ) ;  
 
175
      a.uu = vcat( a.uu , hcat( hcat(bloc(2,0).uu,bloc(2,1).uu) , bloc(2,2).uu ) ) ;  
 
176
      a.ub = vcat( a.ub , hcat( hcat(bloc(2,0).ub,bloc(2,1).ub) , bloc(2,2).ub ) ) ;  
 
177
      a.bu = vcat( a.bu , hcat( hcat(bloc(2,0).bu,bloc(2,1).bu) , bloc(2,2).bu ) ) ;  
 
178
      a.bb = vcat( a.bb , hcat( hcat(bloc(2,0).bb,bloc(2,1).bb) , bloc(2,2).bb ) ) ;  
 
179
      break ;
 
180
    case 6:
 
181
      // tenseur(3D) -> vecteur(3D)
 
182
      break ;
 
183
    default:
 
184
      error_macro("vecteur (" << _ncolbloc << "D) -> tenseur (2D)");
 
185
    }      
 
186
    break ;
 
187
    
 
188
  case 6:
 
189
    trace_macro("3D not yet");
 
190
    switch(_ncolbloc) {
 
191
    case 3: break ;// vecteur(3D) -> tenseur(3D)
 
192
    case 6: break ;// tenseur(3D) -> tenseur(3D)
 
193
    default:
 
194
      error_macro("vecteur (" << _ncolbloc << "D) -> tenseur (3D)");
 
195
    }
 
196
    break ;
 
197
 
 
198
  default:
 
199
    error_macro("vector (" << _ncolbloc << "D) -> vector (" << _nrowbloc << "D)");
 
200
  } // fin du switch(_nrowbloc)
 
201
  return (*this) ;
 
202
}
 
203
form_manip&
 
204
form_manip::operator << (form_manip& (*pf)(form_manip&))
 
205
{
 
206
    return (*pf)(*this);
 
207
}  
 
208
form_manip&
 
209
verbose (form_manip &fm) 
 
210
{
 
211
    fm._verbose = true;
 
212
    return fm;
 
213
}
 
214
form_manip&
 
215
noverbose (form_manip &fm) 
 
216
{
 
217
    fm._verbose = false;
 
218
    return fm;
 
219
}
 
220
size::size (size_type nrowbloc1, size_type ncolbloc1)
 
221
  : _nrowbloc(nrowbloc1),
 
222
    _ncolbloc(ncolbloc1)
 
223
{
 
224
}
 
225
form_manip&
 
226
operator << (form_manip &fm, const size& s)
 
227
{
 
228
    fm._nrowbloc = s._nrowbloc;
 
229
    fm._ncolbloc = s._ncolbloc;
 
230
    return fm ;
 
231
}