AFEPack
BoundaryCondition.h
浏览该文件的文档。
00001 
00011 #ifndef __BoundaryCondition_h__
00012 #define __BoundaryCondition_h__
00013 
00014 #include <cstdarg>
00015 #include <vector>
00016 #include <set>
00017 #include <map>
00018 
00019 #include <lac/sparse_matrix.h>
00020 #include <lac/vector.h>
00021 
00022 #include "Miscellaneous.h"
00023 
00037 class BCondition {
00038  public:
00040   static const int DIRICHLET;
00041   static const int NEUMANN;
00042   static const int ROBIN;
00043 
00044  public:
00045   /*@\{*/
00047   BCondition() {}
00048   BCondition(int type) : _type(type) {}
00049   BCondition(const BCondition& b) :
00050     _type(b._type), _bmark(b._bmark) {}
00051   virtual ~BCondition() {}
00052   /*@\}*/
00053 
00054  public:
00056   BCondition& operator=(const BCondition& b) {
00057     _type = b._type;
00058     _bmark = b._bmark;
00059 
00060     return *this;
00061   }
00062 
00063   /*@\{*/
00065   int type() const {return _type;}
00066   int& type() {return _type;}
00067   /*@\}*/
00068 
00069   /*@\{*/
00071   const std::set<int, std::less<int> >& bound_mark() const {return _bmark;}
00072   std::set<int, std::less<int> >& bound_mark() {return _bmark;}
00073   /*@\}*/
00074 
00075   /*@\{*/
00077   void add_mark(int bm) {add_one_mark(bm);}
00078   void add_mark(const std::vector<int>& bm) {
00079     u_int n = bm.size();
00080     for (u_int i = 0;i < n;++ i) {
00081       add_one_mark(bm[i]);
00082     }
00083   }
00084   template <class V>
00085     void add_mark(u_int n, const V& bm) {
00086     for (u_int i = 0;i < n;i ++) {
00087       add_one_mark(bm[i]);
00088     }
00089   }
00090   void add_mark(u_int n, int bm0, int bm1, ...) {
00091     add_one_mark(bm0);
00092     add_one_mark(bm1);
00093 
00094     va_list ap;
00095     va_start(ap, bm1);
00096     for (u_int i = 2;i < n;i ++) {
00097       add_one_mark(va_arg(ap, int));
00098     }
00099     va_end(ap);
00100   }
00101   /*@\}*/
00102 
00103  private:
00104   void add_one_mark(int bm) {
00105     if (bm == 0) {
00106       std::cerr << "棺ϱʶ 0 ʾڲܼӸ߽"
00107                 << std::endl;
00108       return;
00109     }
00110     _bmark.insert(bm);
00111   }
00112 
00113  public:
00115   virtual void value(const void * p, void * v) const {}
00117   virtual void gradient(const void * p, void * g) const {}
00118 
00119  private:
00120   int _type; 
00121   std::set<int, std::less<int> > _bmark; 
00122 };
00123 
00124 
00125 
00130 template <class P, class V, class G = std::vector<V> >
00131 class BCFunction : public BCondition {
00132   public:
00133   typedef void (*val_fun_t)(const P&, V&);
00134   typedef void (*grad_fun_t)(const P&, G&);
00135 
00136   private:
00137   static void _default_val(const P& p, V& v) {v = V();}
00138   static void _default_grad(const P& p, G& g) {g = G();}
00139 
00140   val_fun_t _val_fun; 
00141   grad_fun_t _grad_fun; 
00142 
00143   public:
00144   BCFunction(val_fun_t vf = &_default_val, 
00145              grad_fun_t gf = &_default_grad) :
00146   _val_fun(vf), 
00147   _grad_fun(gf) {}
00148 
00149   BCFunction(int type, 
00150              val_fun_t vf = &_default_val, 
00151              grad_fun_t gf = &_default_grad) : 
00152   BCondition(type), 
00153   _val_fun(vf), 
00154   _grad_fun(gf) {}
00155 
00156   virtual ~BCFunction() {}
00157 
00158   public:
00159   /*@\{*/
00161   const val_fun_t& value_fun_ptr() const {
00162     return _val_fun;
00163   }
00164   val_fun_t& value_fun_ptr() {
00165     return _val_fun;
00166   }
00167 
00168   const grad_fun_t& gradient_fun_ptr() const {
00169     return _grad_fun;
00170   }
00171   grad_fun_t& gradient_fun_ptr() {
00172     return _grad_fun;
00173   }
00174   /*@\}*/
00175 
00176   public:
00177   /*@\{*/
00179   virtual void value(const void * p, void * v) const {
00180     (*_val_fun)(*(const P *)p, *(V *)v);
00181   }
00182   virtual void gradient(const void * p, void * g) const {
00183     (*_grad_fun)(*(const P *)p, *(G *)g);
00184   }
00185   /*@\}*/
00186 };
00187 
00188 
00216 class BCAdmin {
00217  public:
00218   typedef BCondition * bc_ptr_t;
00219 
00220  private:
00221   std::map<int, bc_ptr_t, std::less<int> > _map;
00222 
00223  public:
00228   template <
00229     class SP,    
00230     class MAT,   
00231     class XVEC,  
00232     class BVEC   
00233     >
00234     void apply(const SP& sp,
00235                MAT& A, 
00236                XVEC& u, 
00237                BVEC& f, 
00238                bool preserve_symmetry = true) 
00239     {
00240       u_int n_dof = sp.n_dof();
00241       const SparsityPattern& spA = A.get_sparsity_pattern();
00242       const std::size_t * rowstart = spA.get_rowstart_indices();
00243       const u_int * colnum = spA.get_column_numbers();
00244       for (u_int i = 0;i < n_dof;++ i) {
00245         int bm = sp.dofInfo(i).boundary_mark;
00246         if (bm == 0) continue; 
00247         const bc_ptr_t bc = this->find(bm);
00249         if (bc == NULL) continue;
00250         if (bc->type() != BCondition::DIRICHLET) continue;
00251         bc->value((const void *)(&sp.dofInfo(i).interp_point), 
00252                   (void *)(&u(i)));
00253 
00254         f(i) = A.diag_element(i)*u(i);
00255         for (u_int j = rowstart[i] + 1;j < rowstart[i + 1];++ j) {
00256           A.global_entry(j) -= A.global_entry(j);
00257         }
00258         if (preserve_symmetry) {
00259           for (u_int j = rowstart[i] + 1;j < rowstart[i + 1];++ j) {
00260             u_int k = colnum[j];
00261             const u_int * p = std::find(&colnum[rowstart[k] + 1],
00262                                         &colnum[rowstart[k + 1]], i);
00263             if (p != &colnum[rowstart[k+1]]) {
00264               u_int l = p - &colnum[rowstart[0]];
00265               f(k) -= A.global_entry(l)*f(i)/A.diag_element(i);
00266               A.global_entry(l) -= A.global_entry(l);
00267             }
00268           }
00269         }
00270       }
00271     }
00272 
00277   template <
00278     class SP,  
00279     class VEC  
00280     >
00281     void clear_entry(const SP& sp,
00282                      VEC& f)
00283     {
00284       u_int n_dof = sp.n_dof();
00285       for (u_int i = 0;i < n_dof;++ i) {
00286         int bm = sp.dof_info(i).bound_mark();
00287         if (bm == 0) continue; 
00288         const bc_ptr_t bc = find(bm);
00289         if (bc != NULL) f(i) -= f(i);
00290       }
00291     }
00292 
00299   template <class BC> 
00300     void add(BC& b) 
00301     {
00302       const std::set<int, std::less<int> >& bm = b.bound_mark();
00303       std::set<int, std::less<int> >::const_iterator 
00304         the_bm = bm.begin(), end_bm = bm.end();
00305       for (;the_bm != end_bm;++ the_bm) {
00306         _map[*the_bm] = &b;
00307       }
00308     }
00309 
00315   bc_ptr_t find(int bm) const 
00316     {
00317       std::map<int, bc_ptr_t, std::less<int> >::const_iterator
00318         the_ptr = _map.find(bm);
00319       if (the_ptr != _map.end()) {
00320         return the_ptr->second;
00321       } else {
00322         return NULL;
00323       }
00324     }
00325 };
00326 
00327 #endif // __BoundaryCondition_h__
00328