00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #ifndef EPOT_PROBLEM_HPP
00044 #define EPOT_PROBLEM_HPP 1
00045
00046
00047 #include <iostream>
00048 #include <stdint.h>
00049 #include "problem.hpp"
00050 #include "solver.hpp"
00051 #include "scalarfield.hpp"
00052 #include "geometry.hpp"
00053 #include "vec3d.hpp"
00054
00055
00068 enum plasma_mode_e {PLASMA_NONE = 0, PLASMA_PEXP_INITIAL, PLASMA_PEXP,
00069 PLASMA_NSIMP_INITIAL, PLASMA_NSIMP};
00070
00071 #define PLASMA_INITIAL PLASMA_PEXP_INITIAL
00072
00142 class EpotProblem : public Problem {
00143
00153 class Node2DoF {
00154 Int3D _size;
00155 int32_t *_n2d;
00157 public:
00158
00159 Node2DoF() : _size(0), _n2d(0) {}
00160 Node2DoF( Int3D size ) : _size(size) {
00161 _n2d = new int32_t[_size[0]*_size[1]*_size[2]];
00162 }
00163 ~Node2DoF() { delete _n2d; }
00164
00165 void resize( Int3D size ) {
00166 _size = size;
00167 if( _n2d )
00168 delete _n2d;
00169 _n2d = new int32_t[_size[0]*_size[1]*_size[2]];
00170 }
00171
00172 int32_t &operator()( int i )
00173 { return( _n2d[i] ); }
00174 int32_t &operator()( int i, int j )
00175 { return( _n2d[i+j*_size[0]] ); }
00176 int32_t &operator()( int i, int j, int k )
00177 { return( _n2d[i+j*_size[0]+k*_size[0]*_size[1]] ); }
00178
00179 const int32_t &operator()( int i ) const
00180 { return( _n2d[i] ); }
00181 const int32_t &operator()( int i, int j ) const
00182 { return( _n2d[i+j*_size[0]] ); }
00183 const int32_t &operator()( int i, int j, int k ) const
00184 { return( _n2d[i+j*_size[0]+k*_size[0]*_size[1]] ); }
00185
00186 void debug_print( void ) const;
00187 };
00188
00189 int32_t _nodecount;
00190 int32_t _dof;
00191 Node2DoF _n2d;
00192 CRowMatrix *_fd_mat;
00193 Vector *_fd_vec;
00195 const Geometry *_g;
00196 mutable CRowMatrix *_fd_mat2;
00197 mutable Vector *_fd_vec2;
00198 mutable Vector *_fd_vec3;
00200 int32_t _neumann_order;
00201 bool _smooth_solid;
00203 plasma_mode_e _plasma;
00205 double _rhoe;
00206 double _Te;
00207 double _Up;
00209 std::vector<double> _rhoi;
00211 std::vector<double> _Ei;
00214 double _force_pot;
00216 bool (*_force_pot_func)(double,double,double);
00217 bool (*_init_plasma_func)(double,double,double);
00219 Solver *_solver;
00222 void set_link( CRowMatrix &A, Vector &B,
00223 int32_t a, int32_t b, double val );
00224
00225 void add_initial_plasma( int32_t i, int32_t j, int32_t k,
00226 CRowMatrix &A, Vector &B, Node2DoF &n2d );
00227
00228 void add_forced_pot( int32_t i, int32_t j, int32_t k,
00229 CRowMatrix &A, Vector &B, Node2DoF &n2d );
00230
00231 void add_vacuum_node( int32_t i, int32_t j, int32_t k,
00232 CRowMatrix &A, Vector &B, Node2DoF &n2d );
00233
00234 void add_neumann_node( signed char a, int32_t i, int32_t j, int32_t k,
00235 CRowMatrix &A, Vector &B, Node2DoF &n2d );
00236
00237 void add_solid_edge_node( signed char a, int32_t i, int32_t j, int32_t k,
00238 CRowMatrix &A, Vector &B, Node2DoF &n2d );
00239
00240 void clear_problem( void );
00241
00242 public:
00243
00244
00245
00246
00247
00250 EpotProblem();
00251
00254 EpotProblem( std::istream &s );
00255
00258 ~EpotProblem();
00259
00260
00261
00262
00263
00268 void set_neumann_order( int32_t order );
00269
00274 void enable_smooth_solids( bool enable );
00275
00283 void set_forced_potential_volume( double force_pot,
00284 bool (*force_pot_func)(double,double,double) );
00285
00290 void set_initial_plasma( double Up,
00291 bool (*plasma_func)(double,double,double) );
00292
00295 void set_pexp_plasma( double rhoe, double Te, double Up );
00296
00302 void set_nsimp_initial_plasma( bool (*plasma_func)(double,double,double) );
00303
00317 void set_nsimp_plasma( double rhop, double Ep,
00318 std::vector<double> rhoi, std::vector<double> Ei );
00319
00324 void construct( const Geometry &g );
00325
00328 void set_solver( Solver &s );
00329
00337 void solve( ScalarField &epot, const ScalarField &scharge ) const;
00338
00339
00340
00341
00342
00350 void get_vecmat( const Matrix **A, const Vector **B ) const;
00351
00359 void get_resjac( const Matrix **J, const Vector **R, const Vector &X ) const;
00360
00363 bool linear( void ) const;
00364
00365
00366
00367
00368
00371 int get_dof( void ) const { return( _dof ); }
00372
00375 void debug_print( void ) const;
00376
00379 void save( std::ostream &s ) const;
00380 };
00381
00382
00383 #endif
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402